Line data Source code
1 : /* Copyright (C) 2000 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 : /********************************************************************/
16 : /** **/
17 : /** GENERIC OPERATIONS **/
18 : /** (first part) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_mod
25 :
26 : /* assume z[1] was created last */
27 : #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
28 : z = gerepileupto((pari_sp)(z+3), gel(z,1));
29 :
30 : /* assume z[1] was created last */
31 : #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
32 : z = gerepileupto((pari_sp)(z+3), gel(z,1));\
33 : else\
34 : gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
35 :
36 : static void
37 105 : warn_coercion(GEN x, GEN y, GEN z)
38 : {
39 105 : if (DEBUGLEVEL)
40 56 : pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
41 105 : }
42 :
43 : static long
44 28 : kro_quad(GEN x, GEN y)
45 28 : { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
46 :
47 : /* is -1 not a square in Zp, assume p prime */
48 : INLINE int
49 42 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
50 :
51 : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
52 : static GEN mulpp(GEN x, GEN y);
53 : static GEN divpp(GEN x, GEN y);
54 : /* Argument codes for inline routines
55 : * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
56 : * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
57 : * T: some type (to be converted to PADIC)
58 : */
59 : static GEN
60 303926510 : addRc(GEN x, GEN y) {
61 303926510 : GEN z = cgetg(3,t_COMPLEX);
62 303907846 : gel(z,1) = gadd(x,gel(y,1));
63 303915515 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 351212885 : mulRc(GEN x, GEN y) {
67 351212885 : GEN z = cgetg(3,t_COMPLEX);
68 351151106 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
69 351560266 : gel(z,2) = gmul(x,gel(y,2)); return z;
70 : }
71 : /* for INTMODs: can't simplify when Re(x) = gen_0 */
72 : static GEN
73 49 : mulRc_direct(GEN x, GEN y) {
74 49 : GEN z = cgetg(3,t_COMPLEX);
75 49 : gel(z,1) = gmul(x,gel(y,1));
76 49 : gel(z,2) = gmul(x,gel(y,2)); return z;
77 : }
78 : static GEN
79 2358035 : divRc(GEN x, GEN y) {
80 2358035 : GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
81 2358027 : GEN z = cgetg(3,t_COMPLEX);
82 2358030 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
83 2358025 : gel(z,2) = gmul(mt, gel(y,2));
84 2358029 : return z;
85 : }
86 : static GEN
87 23128116 : divcR(GEN x, GEN y) {
88 23128116 : GEN z = cgetg(3,t_COMPLEX);
89 23128006 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
90 23127620 : gel(z,2) = gdiv(gel(x,2), y); return z;
91 : }
92 : static GEN
93 1274 : addRq(GEN x, GEN y) {
94 1274 : GEN z = cgetg(4,t_QUAD);
95 1274 : gel(z,1) = ZX_copy(gel(y,1));
96 1274 : gel(z,2) = gadd(x, gel(y,2));
97 1274 : gel(z,3) = gcopy(gel(y,3)); return z;
98 : }
99 : static GEN
100 2310 : mulRq(GEN x, GEN y) {
101 2310 : GEN z = cgetg(4,t_QUAD);
102 2310 : gel(z,1) = ZX_copy(gel(y,1));
103 2310 : gel(z,2) = gmul(x,gel(y,2));
104 2310 : gel(z,3) = gmul(x,gel(y,3)); return z;
105 : }
106 : static GEN
107 77 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
108 77 : long i = gexpo(x) - gexpo(y);
109 77 : if (i > 0) prec += nbits2extraprec(i);
110 77 : return gerepileupto(av, gadd(y, quadtofp(x, prec)));
111 : }
112 : static GEN
113 37585225 : mulrfrac(GEN x, GEN y)
114 : {
115 : pari_sp av;
116 37585225 : GEN z, a = gel(y,1), b = gel(y,2);
117 37585225 : if (is_pm1(a)) /* frequent special case */
118 : {
119 12720567 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
120 12720574 : return z;
121 : }
122 24864419 : av = avma;
123 24864419 : return gerepileuptoleaf(av, divri(mulri(x,a), b));
124 : }
125 : static GEN
126 28 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
127 28 : return gerepileupto(av, gmul(y, quadtofp(x, prec)));
128 : }
129 : static GEN
130 63 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
131 63 : return gerepileupto(av, gdiv(quadtofp(x,prec), y));
132 : }
133 : static GEN
134 42 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
135 42 : return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
136 : }
137 : /* y PADIC, x + y by converting x to padic */
138 : static GEN
139 7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
140 :
141 7 : if (!valp(y)) z = cvtop2(x,y);
142 : else {
143 7 : long l = signe(padic_u(y))? valp(y) + precp(y): valp(y);
144 7 : z = cvtop(x, padic_p(y), l);
145 : }
146 7 : return gerepileupto(av, addsub_pp(z, y, addii));
147 : }
148 : /* y PADIC, x * y by converting x to padic */
149 : static GEN
150 4982871 : mulTp(GEN x, GEN y) { pari_sp av = avma;
151 4982871 : return gerepileupto(av, mulpp(cvtop2(x,y), y));
152 : }
153 : /* y PADIC, non zero x / y by converting x to padic */
154 : static GEN
155 3660 : divTp(GEN x, GEN y) { pari_sp av = avma;
156 3660 : return gerepileupto(av, divpp(cvtop2(x,y), y));
157 : }
158 : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
159 : * converted to O(p^e) and division by 0 */
160 : static GEN
161 1204910 : divpT(GEN x, GEN y) { pari_sp av = avma;
162 1204910 : return gerepileupto(av, divpp(x, cvtop2(y,x)));
163 : }
164 :
165 : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
166 : * clean memory from z on */
167 : static GEN
168 3550198 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
169 3550198 : if (lgefint(X) == 3) {
170 3227073 : ulong u = Fl_add(itou(x),itou(y), X[2]);
171 3227073 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
172 : }
173 : else {
174 323125 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
175 323125 : gel(z,2) = gc_INT((pari_sp)z, u);
176 : }
177 3550198 : gel(z,1) = icopy(X); return z;
178 : }
179 : static GEN
180 1053721 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 1053721 : if (lgefint(X) == 3) {
182 679687 : ulong u = Fl_sub(itou(x),itou(y), X[2]);
183 679687 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
184 : }
185 : else {
186 374034 : GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
187 374034 : gel(z,2) = gc_INT((pari_sp)z, u);
188 : }
189 1053721 : gel(z,1) = icopy(X); return z;
190 : }
191 : /* cf add_intmod_same */
192 : static GEN
193 2835666 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
194 2835666 : if (lgefint(X) == 3) {
195 2694334 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
196 2694334 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
197 : }
198 : else
199 141332 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x,y), X) );
200 2835669 : gel(z,1) = icopy(X); return z;
201 : }
202 : /* cf add_intmod_same */
203 : static GEN
204 341973 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
205 : {
206 341973 : if (lgefint(X) == 3) {
207 319431 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
208 319368 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 22542 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
212 341910 : gel(z,1) = icopy(X); return z;
213 : }
214 :
215 : /*******************************************************************/
216 : /* */
217 : /* REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC) */
218 : /* */
219 : /* (static routines are not memory clean, but OK for gerepileupto) */
220 : /*******************************************************************/
221 : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
222 : * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
223 : static GEN
224 9858014 : rfrac_denom_mul_scal(GEN d, GEN y)
225 : {
226 9858014 : GEN D = RgX_Rg_mul(d, y);
227 9858014 : if (lg(D) != lg(d))
228 : { /* try to generate a meaningful diagnostic */
229 0 : D = gdiv(leading_coeff(d), y); /* should fail */
230 0 : pari_err_INV("gred_rfrac", y); /* better than nothing */
231 : }
232 9858014 : return D;
233 : }
234 :
235 : static int
236 57698360 : Leading_is_neg(GEN x)
237 : {
238 122017057 : while (typ(x) == t_POL) x = leading_coeff(x);
239 57698360 : return is_real_t(typ(x))? gsigne(x) < 0: 0;
240 : }
241 :
242 : static int
243 154191271 : transtype(GEN x) { return x != gen_1 && typ(x) != t_PADIC; }
244 :
245 : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
246 : GEN
247 57702596 : gred_rfrac_simple(GEN n, GEN d)
248 : {
249 : GEN _1n, _1d, c, cn, cd, z;
250 57702596 : long dd = degpol(d);
251 :
252 57702596 : if (dd <= 0)
253 : {
254 4236 : if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
255 4236 : n = gdiv(n, gel(d,2));
256 4236 : if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
257 4236 : return n;
258 : }
259 57698360 : if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
260 57698360 : _1n = Rg_get_1(n);
261 57698360 : _1d = Rg_get_1(d);
262 57698360 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
263 57698360 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
264 57698353 : cd = content(d);
265 59564964 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
266 57698353 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
267 57698353 : if (!gequal1(cd)) {
268 6568219 : d = RgX_Rg_div(d,cd);
269 6568219 : if (!gequal1(cn))
270 : {
271 1294266 : if (gequal0(cn)) {
272 49 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
273 0 : n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
274 0 : c = gen_1;
275 : } else {
276 1294217 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
277 1294217 : c = gdiv(cn,cd);
278 : }
279 : }
280 : else
281 5273953 : c = ginv(cd);
282 : } else {
283 51130134 : if (!gequal1(cn))
284 : {
285 3277648 : if (gequal0(cn)) {
286 987 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
287 21 : c = gen_1;
288 : } else {
289 3276661 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
290 3276661 : c = cn;
291 : }
292 : } else {
293 47852486 : GEN y = cgetg(3,t_RFRAC);
294 47852486 : gel(y,1) = gcopy(n);
295 47852486 : gel(y,2) = RgX_copy(d); return y;
296 : }
297 : }
298 :
299 9844852 : if (typ(c) == t_POL)
300 : {
301 901227 : z = c;
302 939657 : do { z = content(z); } while (typ(z) == t_POL);
303 901227 : cd = denom_i(z);
304 901227 : cn = gmul(c, cd);
305 : }
306 : else
307 : {
308 8943625 : cn = numer_i(c);
309 8943625 : cd = denom_i(c);
310 : }
311 9844852 : z = cgetg(3,t_RFRAC);
312 9844852 : gel(z,1) = gmul(n, cn);
313 9844852 : gel(z,2) = d = rfrac_denom_mul_scal(d, cd);
314 : /* can happen: Pol(O(17^-1)) / Pol([Mod(9,23), O(23^-3)]) */
315 9844852 : if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
316 9844852 : return z;
317 : }
318 :
319 : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
320 : static GEN
321 207774 : fix_rfrac(GEN x, long d)
322 : {
323 : GEN z, N, D;
324 207774 : if (!d || typ(x) == t_POL) return x;
325 165872 : z = cgetg(3, t_RFRAC);
326 165872 : N = gel(x,1);
327 165872 : D = gel(x,2);
328 165872 : if (d > 0) {
329 6265 : gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
330 160447 : : monomialcopy(N,d,varn(D));
331 160384 : gel(z, 2) = RgX_copy(D);
332 : } else {
333 5488 : gel(z, 1) = gcopy(N);
334 5488 : gel(z, 2) = RgX_shift(D, -d);
335 : }
336 165872 : return z;
337 : }
338 :
339 : /* assume d != 0 */
340 : static GEN
341 44760592 : gred_rfrac2(GEN n, GEN d)
342 : {
343 : GEN y, z, _1n, _1d;
344 : long v, vd, vn;
345 :
346 44760592 : n = simplify_shallow(n);
347 44760592 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
348 37919688 : d = simplify_shallow(d);
349 37919688 : if (typ(d) != t_POL) return gdiv(n,d);
350 36740944 : vd = varn(d);
351 36740944 : if (typ(n) != t_POL)
352 : {
353 20620789 : if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
354 20619382 : if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
355 0 : pari_err_BUG("gred_rfrac2 [incompatible variables]");
356 : }
357 16120155 : vn = varn(n);
358 16120155 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
359 15982446 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
360 15822673 : _1n = Rg_get_1(n);
361 15822673 : _1d = Rg_get_1(d);
362 15822673 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
363 15822673 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
364 :
365 : /* now n and d are t_POLs in the same variable */
366 15822673 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
367 15822673 : if (!degpol(d))
368 : {
369 12784971 : n = RgX_Rg_div(n,gel(d,2));
370 12784971 : return v? RgX_mulXn(n,v): n;
371 : }
372 :
373 : /* X does not divide gcd(n,d), deg(d) > 0 */
374 3037702 : if (!isinexact(n) && !isinexact(d))
375 : {
376 3037478 : y = RgX_divrem(n, d, &z);
377 3037478 : if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
378 207550 : z = RgX_gcd(d, z);
379 207550 : if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
380 : }
381 207774 : return fix_rfrac(gred_rfrac_simple(n,d), v);
382 : }
383 :
384 : /* x,y t_INT, return x/y in reduced form */
385 : GEN
386 139134232 : Qdivii(GEN x, GEN y)
387 : {
388 139134232 : pari_sp av = avma;
389 : GEN r, q;
390 :
391 139134232 : if (lgefint(y) == 3)
392 : {
393 121149361 : q = Qdiviu(x, y[2]);
394 121147254 : if (signe(y) > 0) return q;
395 11016120 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
396 11016312 : return q;
397 : }
398 17984871 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
399 17986315 : if (equali1(x))
400 : {
401 5182068 : if (!signe(y)) pari_err_INV("gdiv",y);
402 5181963 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
403 : }
404 12804276 : q = dvmdii(x,y,&r);
405 12804183 : if (r == gen_0) return q; /* gen_0 intended */
406 8378723 : r = gcdii(y, r);
407 8378731 : if (lgefint(r) == 3)
408 : {
409 7494616 : ulong t = r[2];
410 7494616 : set_avma(av);
411 7494583 : if (t == 1) q = mkfraccopy(x,y);
412 : else
413 : {
414 2650058 : q = cgetg(3,t_FRAC);
415 2650195 : gel(q,1) = diviuexact(x,t);
416 2650121 : gel(q,2) = diviuexact(y,t);
417 : }
418 : }
419 : else
420 : { /* rare: r and q left on stack for efficiency */
421 884115 : q = cgetg(3,t_FRAC);
422 884123 : gel(q,1) = diviiexact(x,r);
423 884124 : gel(q,2) = diviiexact(y,r);
424 : }
425 8378632 : normalize_frac(q); return q;
426 : }
427 :
428 : /* x t_INT, return x/y in reduced form */
429 : GEN
430 140472984 : Qdiviu(GEN x, ulong y)
431 : {
432 140472984 : pari_sp av = avma;
433 : ulong r, t;
434 : GEN q;
435 :
436 140472984 : if (y == 1) return icopy(x);
437 119690259 : if (!y) pari_err_INV("gdiv",gen_0);
438 119690259 : if (equali1(x)) retmkfrac(gen_1, utoipos(y));
439 112562340 : q = absdiviu_rem(x,y,&r);
440 112558619 : if (!r)
441 : {
442 68859316 : if (signe(x) < 0) togglesign(q);
443 68859326 : return q;
444 : }
445 43699303 : t = ugcd(y, r); set_avma(av);
446 43706060 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
447 17386713 : retmkfrac(diviuexact(x,t), utoipos(y / t));
448 : }
449 :
450 : /* x t_INT, return x/y in reduced form */
451 : GEN
452 1583916 : Qdivis(GEN x, long y)
453 : {
454 1583916 : pari_sp av = avma;
455 : ulong r, t;
456 : long s;
457 : GEN q;
458 :
459 1583916 : if (y > 0) return Qdiviu(x, y);
460 177681 : if (!y) pari_err_INV("gdiv",gen_0);
461 177681 : s = signe(x);
462 177681 : if (!s) return gen_0;
463 132202 : if (y < 0) { y = -y; s = -s; }
464 132202 : if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
465 131922 : if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
466 129759 : q = absdiviu_rem(x,y,&r);
467 129759 : if (!r)
468 : {
469 55356 : if (s < 0) togglesign(q);
470 55356 : return q;
471 : }
472 74403 : t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
473 74403 : if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
474 74403 : gel(q,1) = x; setsigne(x, s);
475 74403 : gel(q,2) = utoipos(y); return q;
476 : }
477 :
478 : /*******************************************************************/
479 : /* */
480 : /* CONJUGATION */
481 : /* */
482 : /*******************************************************************/
483 : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
484 : static GEN
485 18340 : quad_polmod_conj(GEN x, GEN y)
486 : {
487 : GEN z, u, v, a, b;
488 18340 : if (typ(x) != t_POL) return gcopy(x);
489 18340 : if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
490 18340 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
491 18340 : b = gel(y,3); v = gel(x,2);
492 18340 : z = cgetg(4, t_POL); z[1] = x[1];
493 18340 : gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
494 18340 : gel(z,3) = gneg(u); return z;
495 : }
496 : static GEN
497 18340 : quad_polmod_norm(GEN x, GEN y)
498 : {
499 : GEN z, u, v, a, b, c;
500 18340 : if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
501 18340 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
502 18340 : b = gel(y,3); v = gel(x,2);
503 18340 : c = gel(y,2);
504 18340 : z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
505 18340 : if (!gequal1(a)) z = gdiv(z, a);
506 18340 : return gadd(z, gsqr(v));
507 : }
508 :
509 : GEN
510 19577694 : conj_i(GEN x)
511 : {
512 19577694 : switch(typ(x))
513 : {
514 6541247 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
515 6541247 : return x;
516 :
517 12872875 : case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
518 :
519 924 : case t_QUAD:
520 : {
521 924 : GEN y = cgetg(4,t_QUAD);
522 924 : gel(y,1) = gel(x,1);
523 924 : gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
524 924 : : gadd(gel(x,2), gel(x,3));
525 924 : gel(y,3) = gneg(gel(x,3)); return y;
526 : }
527 350 : case t_POL: pari_APPLY_pol_normalized(conj_i(gel(x,i)));
528 31591 : case t_SER: pari_APPLY_ser_normalized(conj_i(gel(x,i)));
529 :
530 152733 : case t_RFRAC:
531 : case t_VEC:
532 : case t_COL:
533 550256 : case t_MAT: pari_APPLY_same(conj_i(gel(x,i)));
534 :
535 0 : case t_POLMOD:
536 : {
537 0 : GEN X = gel(x,1);
538 0 : long d = degpol(X);
539 0 : if (d < 2) return x;
540 0 : if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
541 : }
542 : }
543 0 : pari_err_TYPE("gconj",x);
544 : return NULL; /* LCOV_EXCL_LINE */
545 : }
546 : GEN
547 133221 : gconj(GEN x)
548 133221 : { pari_sp av = avma; return gc_GEN(av, conj_i(x)); }
549 :
550 : GEN
551 84 : conjvec(GEN x,long prec)
552 : {
553 : long lx, s, i;
554 : GEN z;
555 :
556 84 : switch(typ(x))
557 : {
558 0 : case t_INT: case t_INTMOD: case t_FRAC:
559 0 : return mkcolcopy(x);
560 :
561 0 : case t_COMPLEX: case t_QUAD:
562 0 : z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
563 :
564 28 : case t_FFELT:
565 28 : return FF_conjvec(x);
566 :
567 0 : case t_VEC: case t_COL:
568 0 : lx = lg(x); z = cgetg(lx,t_MAT);
569 0 : if (lx == 1) return z;
570 0 : gel(z,1) = conjvec(gel(x,1),prec);
571 0 : s = lgcols(z);
572 0 : for (i=2; i<lx; i++)
573 : {
574 0 : gel(z,i) = conjvec(gel(x,i),prec);
575 0 : if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
576 : }
577 0 : break;
578 :
579 56 : case t_POLMOD: {
580 56 : GEN T = gel(x,1), r;
581 : pari_sp av;
582 :
583 56 : lx = lg(T);
584 56 : if (lx <= 3) return cgetg(1,t_COL);
585 56 : x = gel(x,2);
586 238 : for (i=2; i<lx; i++)
587 : {
588 189 : GEN c = gel(T,i);
589 189 : switch(typ(c)) {
590 7 : case t_INTMOD: {
591 7 : GEN p = gel(c,1);
592 : pari_sp av;
593 7 : if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
594 7 : av = avma;
595 7 : T = RgX_to_FpX(T,p);
596 7 : x = RgX_to_FpX(x, p);
597 7 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
598 7 : z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
599 7 : return gerepileupto(av, z);
600 : }
601 182 : case t_INT:
602 182 : case t_FRAC: break;
603 0 : default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
604 : }
605 : }
606 49 : if (typ(x) != t_POL)
607 : {
608 0 : if (!is_rational_t(typ(x)))
609 0 : pari_err_TYPE("conjvec [not a rational t_POL]",x);
610 0 : retconst_col(lx-3, gcopy(x));
611 : }
612 49 : RgX_check_QX(x,"conjvec");
613 49 : av = avma;
614 49 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
615 49 : r = cleanroots(T,prec);
616 49 : z = cgetg(lx-2,t_COL);
617 182 : for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
618 49 : return gerepileupto(av, z);
619 : }
620 :
621 0 : default:
622 0 : pari_err_TYPE("conjvec",x);
623 : return NULL; /* LCOV_EXCL_LINE */
624 : }
625 0 : return z;
626 : }
627 :
628 : /********************************************************************/
629 : /** **/
630 : /** ADDITION **/
631 : /** **/
632 : /********************************************************************/
633 : static GEN
634 24587927 : mkpadic_mod(GEN u, GEN p, GEN pd, long e, long d)
635 24587927 : { retmkpadic(modii(u, pd), icopy(p), icopy(pd), e, d); }
636 :
637 : /* x, y compatible PADIC, op = add or sub */
638 : static GEN
639 19967403 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
640 : {
641 19967403 : pari_sp av = avma;
642 : long d, e, r, rx, ry;
643 19967403 : GEN u, z, mod, pdx, pdy, ux, uy, p = padic_p(x);
644 : int swap;
645 :
646 19967403 : e = valp(x);
647 19967403 : r = valp(y); d = r-e;
648 19967403 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
649 19967403 : pdx = padic_pd(x); ux = padic_u(x);
650 19967403 : pdy = padic_pd(y); uy = padic_u(y);
651 19967403 : rx = precp(x);
652 19967403 : ry = precp(y);
653 19967403 : if (d) /* v(x) < v(y) */
654 : {
655 10703304 : r = d+ry; z = powiu(p,d);
656 10703472 : if (r < rx) mod = mulii(z,pdy); else { r = rx; mod = pdx; }
657 10703498 : z = mulii(z, uy);
658 10703400 : u = swap? op(z, ux): op(ux, z);
659 : }
660 : else
661 : {
662 : long c;
663 9264099 : if (ry < rx) { r=ry; mod = pdy; } else { r=rx; mod = pdx; }
664 9264099 : u = swap? op(uy, ux): op(ux, uy);
665 9265111 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
666 : {
667 138223 : set_avma(av); return zeropadic(p, e+r);
668 : }
669 9126722 : if (c)
670 : {
671 3351445 : mod = diviiexact(mod, powiu(p,c));
672 3351446 : r -= c;
673 3351446 : e += c;
674 : }
675 : }
676 19830083 : return gerepileupto(av, mkpadic_mod(u, p, mod, e, r));
677 : }
678 : /* Rg_to_Fp(t_FRAC) without GC */
679 : static GEN
680 226838 : Q_to_Fp(GEN x, GEN p)
681 226838 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
682 : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
683 : static GEN
684 4762885 : addQp(GEN x, GEN y)
685 : {
686 4762885 : pari_sp av = avma;
687 4762885 : long d, r, e, vy = valp(y), py = precp(y);
688 4762885 : GEN q, mod, u, p = padic_p(y);
689 :
690 4762885 : e = Q_pvalrem(x, p, &x);
691 4762865 : d = vy - e; r = d + py;
692 4762865 : if (r <= 0) { set_avma(av); return gcopy(y); }
693 4760984 : mod = padic_pd(y);
694 4760984 : u = padic_u(y);
695 4760984 : if (d > 0)
696 : {
697 1376569 : q = powiu(p,d);
698 1376569 : mod = mulii(mod, q);
699 1376567 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
700 1376567 : u = addii(x, mulii(u, q));
701 : }
702 3384415 : else if (d < 0)
703 : {
704 405332 : q = powiu(p,-d);
705 405332 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
706 405332 : u = addii(u, mulii(x, q));
707 405332 : r = py; e = vy;
708 : }
709 : else
710 : {
711 : long c;
712 2979083 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
713 2979083 : u = addii(u, x);
714 2979078 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
715 : {
716 1205 : set_avma(av); return zeropadic(p,e+r);
717 : }
718 2977871 : if (c)
719 : {
720 970644 : mod = diviiexact(mod, powiu(p,c));
721 970645 : r -= c;
722 970645 : e += c;
723 : }
724 : }
725 4759774 : return gerepileupto(av, mkpadic_mod(u, p, mod, e, r));
726 : }
727 :
728 : /* Mod(x,X) + Mod(y,X) */
729 : #define addsub_polmod_same addsub_polmod_scal
730 : /* Mod(x,X) +/- Mod(y,Y) */
731 : static GEN
732 7203 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
733 : {
734 7203 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
735 7203 : GEN z = cgetg(3,t_POLMOD);
736 7203 : long vx = varn(X), vy = varn(Y);
737 7203 : if (vx==vy) {
738 : pari_sp av;
739 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
740 14 : warn_coercion(X,Y,gel(z,1));
741 14 : gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
742 : }
743 7189 : if (varncmp(vx, vy) < 0)
744 7189 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
745 : else
746 0 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
747 7189 : gel(z,2) = op(x, y); return z;
748 : }
749 : /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
750 : static GEN
751 13420724 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
752 13420724 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
753 :
754 : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
755 : static GEN
756 408094 : add_ser_scal(GEN y, GEN x)
757 : {
758 : long i, v, ly, vy;
759 : GEN z;
760 :
761 408094 : if (isrationalzero(x)) return gcopy(y);
762 381837 : ly = lg(y);
763 381837 : v = valser(y);
764 381837 : if (v < 3-ly) return gcopy(y);
765 : /* v + ly >= 3 */
766 381585 : if (v < 0)
767 : {
768 1162 : z = cgetg(ly,t_SER); z[1] = y[1];
769 3276 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
770 1162 : gel(z,i) = gadd(x,gel(y,i)); i++;
771 3157 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
772 1162 : return normalizeser(z);
773 : }
774 380423 : vy = varn(y);
775 380423 : if (v > 0)
776 : {
777 19299 : if (ser_isexactzero(y))
778 9338 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
779 9961 : y -= v; ly += v;
780 9961 : z = cgetg(ly,t_SER);
781 9961 : x = gcopy(x);
782 20559 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
783 : }
784 361124 : else if (ser_isexactzero(y)) return gcopy(y);
785 : else
786 : { /* v = 0, ly >= 3 */
787 361117 : z = cgetg(ly,t_SER);
788 361117 : x = gadd(x, gel(y,2));
789 361117 : i = 3;
790 : }
791 1598979 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
792 371078 : gel(z,2) = x;
793 371078 : z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
794 371078 : return gequal0(x)? normalizeser(z): z;
795 : }
796 : static long
797 7161978 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
798 : /* x,y t_SER in the same variable: x+y */
799 : static GEN
800 3581388 : ser_add(GEN x, GEN y)
801 : {
802 3581388 : long i, lx,ly, n = valser(y) - valser(x);
803 : GEN z;
804 3581388 : if (n < 0) { n = -n; swap(x,y); }
805 : /* valser(x) <= valser(y) */
806 3581388 : lx = _serprec(x);
807 3581388 : if (lx == 2) /* don't lose type information */
808 : {
809 798 : z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
810 798 : setvalser(z, valser(x)); return z;
811 : }
812 3580590 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
813 3580590 : if (n)
814 : {
815 107782 : if (n+2 > lx) return gcopy(x);
816 107082 : z = cgetg(ly,t_SER);
817 819116 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
818 506624 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
819 : } else {
820 3472808 : z = cgetg(ly,t_SER);
821 17589354 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
822 : }
823 3579890 : z[1] = x[1]; return normalizeser(z);
824 : }
825 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
826 : static GEN
827 8811940 : add_rfrac_scal(GEN y, GEN x)
828 : {
829 : pari_sp av;
830 : GEN n;
831 :
832 8811940 : if (isintzero(x)) return gcopy(y); /* frequent special case */
833 5138175 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
834 5138175 : return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
835 : }
836 :
837 : /* x "scalar", ty != t_MAT and nonscalar */
838 : static GEN
839 44614087 : add_scal(GEN y, GEN x, long ty)
840 : {
841 44614087 : switch(ty)
842 : {
843 39701806 : case t_POL: return RgX_Rg_add(y, x);
844 408066 : case t_SER: return add_ser_scal(y, x);
845 4468855 : case t_RFRAC: return add_rfrac_scal(y, x);
846 0 : case t_COL: return RgC_Rg_add(y, x);
847 35360 : case t_VEC:
848 35360 : if (isintzero(x)) return gcopy(y);
849 175 : break;
850 : }
851 175 : pari_err_TYPE2("+",x,y);
852 : return NULL; /* LCOV_EXCL_LINE */
853 : }
854 :
855 : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
856 : static GEN
857 14436958 : setfrac(GEN z, GEN a, GEN b)
858 : {
859 14436958 : gel(z,1) = icopy_avma(a, (pari_sp)z);
860 14436951 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
861 14437033 : set_avma((pari_sp)gel(z,2)); return z;
862 : }
863 : /* z <- a / (b*Q), (Q,a) = 1 */
864 : static GEN
865 13628232 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
866 : {
867 13628232 : GEN q = Qdivii(a, b); /* != 0 */
868 13628303 : if (typ(q) == t_INT)
869 : {
870 1950976 : gel(z,1) = gc_INT((pari_sp)Q, q);
871 1950976 : gel(z,2) = Q; return z;
872 : }
873 11677327 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
874 : }
875 : static GEN
876 41396841 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
877 : {
878 41396841 : GEN x1 = gel(x,1), x2 = gel(x,2);
879 41396841 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
880 41396841 : int s = cmpii(x2, y2);
881 :
882 : /* common denominator: (x1 op y1) / x2 */
883 41396842 : if (!s)
884 : {
885 16469382 : pari_sp av = avma;
886 16469382 : return gerepileupto(av, Qdivii(op(x1, y1), x2));
887 : }
888 24927460 : z = cgetg(3, t_FRAC);
889 24927476 : if (s < 0)
890 : {
891 17263307 : Q = dvmdii(y2, x2, &r);
892 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
893 17263260 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
894 9739958 : delta = gcdii(x2,r);
895 : }
896 : else
897 : {
898 7664169 : Q = dvmdii(x2, y2, &r);
899 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
900 7664175 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
901 1559188 : delta = gcdii(y2,r);
902 : }
903 : /* delta = gcd(x2,y2) */
904 11299171 : if (equali1(delta))
905 : { /* numerator is nonzero */
906 8539463 : gel(z,1) = gc_INT((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
907 8539466 : gel(z,2) = mulii(x2,y2); return z;
908 : }
909 2759704 : x2 = diviiexact(x2,delta);
910 2759705 : y2 = diviiexact(y2,delta);
911 2759709 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
912 2759708 : q = dvmdii(n, delta, &r);
913 2759704 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
914 2506531 : r = gcdii(delta, r);
915 2506535 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
916 2506535 : return setfrac(z, n, mulii(mulii(x2, y2), delta));
917 : }
918 :
919 : /* assume x2, y2 are t_POLs in the same variable */
920 : static GEN
921 3018187 : add_rfrac(GEN x, GEN y)
922 : {
923 3018187 : pari_sp av = avma;
924 3018187 : GEN x1 = gel(x,1), x2 = gel(x,2);
925 3018187 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
926 :
927 3018187 : delta = RgX_gcd(x2,y2);
928 3018187 : if (!degpol(delta))
929 : {
930 644 : n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
931 644 : d = RgX_mul(x2, y2);
932 644 : return gerepileupto(av, gred_rfrac_simple(n, d));
933 : }
934 3017543 : x2 = RgX_div(x2,delta);
935 3017543 : y2 = RgX_div(y2,delta);
936 3017543 : n = gadd(gmul(x1,y2), gmul(y1,x2));
937 3017543 : if (!signe(n))
938 : {
939 721418 : n = simplify_shallow(n);
940 721418 : if (isexactzero(n))
941 : {
942 721411 : if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
943 7 : return gerepileupto(av, scalarpol(n, varn(x2)));
944 : }
945 7 : return gc_GEN(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
946 : }
947 2296125 : if (degpol(n) == 0)
948 1150575 : return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
949 1145550 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
950 1145550 : if (isexactzero(r))
951 : {
952 : GEN z;
953 227982 : d = RgX_mul(x2, y2);
954 : /* "constant" denominator ? */
955 227982 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
956 227982 : return gerepileupto(av, z);
957 : }
958 917568 : r = RgX_gcd(delta, r);
959 917568 : if (degpol(r))
960 : {
961 160642 : n = RgX_div(n, r);
962 160642 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
963 : }
964 : else
965 756926 : d = RgX_mul(gel(x,2), y2);
966 917568 : return gerepileupto(av, gred_rfrac_simple(n, d));
967 : }
968 :
969 : GEN
970 6072906686 : gadd(GEN x, GEN y)
971 : {
972 6072906686 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
973 : pari_sp av, tetpil;
974 : GEN z, p1;
975 :
976 6072906686 : if (tx == ty) switch(tx) /* shortcut to generic case */
977 : {
978 2995246665 : case t_INT: return addii(x,y);
979 1638111677 : case t_REAL: return addrr(x,y);
980 1381493 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
981 1381493 : z = cgetg(3,t_INTMOD);
982 1381493 : if (X==Y || equalii(X,Y))
983 1381479 : return add_intmod_same(z, X, gel(x,2), gel(y,2));
984 14 : gel(z,1) = gcdii(X,Y);
985 14 : warn_coercion(X,Y,gel(z,1));
986 14 : av = avma; p1 = addii(gel(x,2),gel(y,2));
987 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
988 : }
989 36325702 : case t_FRAC: return addsub_frac(x,y,addii);
990 354463571 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
991 353985835 : gel(z,2) = gadd(gel(x,2),gel(y,2));
992 354840283 : if (isintzero(gel(z,2)))
993 : {
994 521178 : set_avma((pari_sp)(z+3));
995 521178 : return gadd(gel(x,1),gel(y,1));
996 : }
997 354143152 : gel(z,1) = gadd(gel(x,1),gel(y,1));
998 354267742 : return z;
999 14529468 : case t_PADIC:
1000 14529468 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1001 14529135 : return addsub_pp(x,y, addii);
1002 672 : case t_QUAD: z = cgetg(4,t_QUAD);
1003 672 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1004 672 : gel(z,1) = ZX_copy(gel(x,1));
1005 672 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1006 672 : gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
1007 10156341 : case t_POLMOD:
1008 10156341 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1009 10149173 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
1010 7168 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
1011 13133545 : case t_FFELT: return FF_add(x,y);
1012 139552789 : case t_POL:
1013 139552789 : vx = varn(x);
1014 139552789 : vy = varn(y);
1015 139552789 : if (vx != vy) {
1016 4856925 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1017 1798442 : else return RgX_Rg_add(y, x);
1018 : }
1019 134695864 : return RgX_add(x, y);
1020 3576033 : case t_SER:
1021 3576033 : vx = varn(x);
1022 3576033 : vy = varn(y);
1023 3576033 : if (vx != vy) {
1024 28 : if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1025 21 : else return add_ser_scal(y, x);
1026 : }
1027 3576005 : return ser_add(x, y);
1028 4310088 : case t_RFRAC:
1029 4310088 : vx = varn(gel(x,2));
1030 4310088 : vy = varn(gel(y,2));
1031 4310088 : if (vx != vy) {
1032 1291901 : if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1033 538397 : else return add_rfrac_scal(y, x);
1034 : }
1035 3018187 : return add_rfrac(x,y);
1036 2060826 : case t_VEC:
1037 2060826 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1038 2060826 : return RgV_add(x,y);
1039 1103487 : case t_COL:
1040 1103487 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1041 1103487 : return RgC_add(x,y);
1042 671412 : case t_MAT:
1043 671412 : lx = lg(x);
1044 671412 : if (lg(y) != lx) pari_err_OP("+",x,y);
1045 671412 : if (lx == 1) return cgetg(1, t_MAT);
1046 671412 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1047 671405 : return RgM_add(x,y);
1048 :
1049 0 : default: pari_err_TYPE2("+",x,y);
1050 : }
1051 : /* tx != ty */
1052 862500433 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1053 :
1054 862500433 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1055 : {
1056 746978218 : case t_INT:
1057 : switch(ty)
1058 : {
1059 436728862 : case t_REAL: return addir(x,y);
1060 2139125 : case t_INTMOD:
1061 2139125 : z = cgetg(3, t_INTMOD);
1062 2139125 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1063 51133059 : case t_FRAC: z = cgetg(3,t_FRAC);
1064 51133055 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1065 51132994 : gel(z,2) = icopy(gel(y,2)); return z;
1066 250055945 : case t_COMPLEX: return addRc(x, y);
1067 5597426 : case t_PADIC:
1068 5597426 : if (!signe(x)) return gcopy(y);
1069 4533148 : return addQp(x,y);
1070 1113 : case t_QUAD: return addRq(x, y);
1071 1365100 : case t_FFELT: return FF_Z_add(y,x);
1072 : }
1073 :
1074 : case t_REAL:
1075 57322550 : switch(ty)
1076 : {
1077 13734470 : case t_FRAC:
1078 13734470 : if (!signe(gel(y,1))) return rcopy(x);
1079 13734470 : if (!signe(x))
1080 : {
1081 12080 : lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1082 12080 : return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1083 : }
1084 13722390 : av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
1085 13720984 : return gerepile(av,tetpil,divri(z,gel(y,2)));
1086 43588010 : case t_COMPLEX: return addRc(x, y);
1087 70 : case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, realprec(x));
1088 :
1089 0 : default: pari_err_TYPE2("+",x,y);
1090 : }
1091 :
1092 17640 : case t_INTMOD:
1093 : switch(ty)
1094 : {
1095 17507 : case t_FRAC: { GEN X = gel(x,1);
1096 17507 : z = cgetg(3, t_INTMOD);
1097 17507 : p1 = Fp_div(gel(y,1), gel(y,2), X);
1098 17507 : return add_intmod_same(z, X, p1, gel(x,2));
1099 : }
1100 14 : case t_FFELT:
1101 14 : if (!equalii(gel(x,1),FF_p_i(y)))
1102 0 : pari_err_OP("+",x,y);
1103 14 : return FF_Z_add(y,gel(x,2));
1104 91 : case t_COMPLEX: return addRc(x, y);
1105 0 : case t_PADIC: { GEN X = gel(x,1);
1106 0 : z = cgetg(3, t_INTMOD);
1107 0 : return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1108 : }
1109 28 : case t_QUAD: return addRq(x, y);
1110 : }
1111 :
1112 : case t_FRAC:
1113 : switch (ty)
1114 : {
1115 10299659 : case t_COMPLEX: return addRc(x, y);
1116 229736 : case t_PADIC:
1117 229736 : if (!signe(gel(x,1))) return gcopy(y);
1118 229736 : return addQp(x,y);
1119 133 : case t_QUAD: return addRq(x, y);
1120 1337 : case t_FFELT: return FF_Q_add(y, x);
1121 : }
1122 :
1123 : case t_FFELT:
1124 0 : pari_err_TYPE2("+",x,y);
1125 :
1126 35 : case t_COMPLEX:
1127 : switch(ty)
1128 : {
1129 28 : case t_PADIC:
1130 28 : return Zp_nosquare_m1(padic_p(y))? addRc(y, x): addTp(x, y);
1131 7 : case t_QUAD:
1132 7 : lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1133 7 : return gequal0(y)? gcopy(x): addqf(y, x, lx);
1134 : }
1135 :
1136 : case t_PADIC: /* ty == t_QUAD */
1137 0 : return (kro_quad(y, padic_p(x)) == -1)? addRq(x, y): addTp(y, x);
1138 : }
1139 : /* tx < ty, !is_const_t(y) */
1140 50285405 : switch(ty)
1141 : {
1142 7693 : case t_MAT:
1143 7693 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1144 7693 : if (isrationalzero(x)) return gcopy(y);
1145 7616 : return RgM_Rg_add(y, x);
1146 185394 : case t_COL:
1147 185394 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1148 185394 : return RgC_Rg_add(y, x);
1149 2391602 : case t_POLMOD: /* is_const_t(tx) in this case */
1150 2391602 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1151 : }
1152 47700716 : if (is_scalar_t(tx)) {
1153 44620873 : if (tx == t_POLMOD)
1154 : {
1155 109851 : vx = varn(gel(x,1));
1156 109851 : vy = gvar(y);
1157 109851 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1158 : else
1159 80850 : if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1160 29540 : return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1161 : }
1162 44511022 : return add_scal(y, x, ty);
1163 : }
1164 : /* x and y are not scalars, ty != t_MAT */
1165 3079845 : vx = gvar(x);
1166 3079851 : vy = gvar(y);
1167 3079851 : if (vx != vy) { /* x or y is treated as a scalar */
1168 22759 : if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1169 32417 : return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1170 32417 : : add_scal(y, x, ty);
1171 : }
1172 : /* vx = vy */
1173 3057092 : switch(tx)
1174 : {
1175 3056595 : case t_POL:
1176 : switch (ty)
1177 : {
1178 5411 : case t_SER:
1179 5411 : if (lg(x) == 2) return gcopy(y);
1180 5390 : i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1181 5390 : i = lg(y) + valser(y) - i;
1182 5390 : if (i < 3) return gcopy(y);
1183 5383 : p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1184 5383 : settyp(p1, t_VECSMALL); /* p1 left on stack */
1185 5383 : return y;
1186 :
1187 3051184 : case t_RFRAC: return add_rfrac_scal(y, x);
1188 : }
1189 0 : break;
1190 :
1191 497 : case t_SER:
1192 497 : if (ty == t_RFRAC)
1193 : {
1194 : long vn, vd;
1195 497 : av = avma;
1196 497 : vn = gval(gel(y,1), vy);
1197 497 : vd = RgX_valrem_inexact(gel(y,2), NULL);
1198 497 : if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
1199 :
1200 497 : l = lg(x) + valser(x) - (vn - vd);
1201 497 : if (l < 3) { set_avma(av); return gcopy(x); }
1202 497 : return gerepileupto(av, gadd(x, rfrac_to_ser_i(y, l)));
1203 : }
1204 0 : break;
1205 : }
1206 0 : pari_err_TYPE2("+",x,y);
1207 : return NULL; /* LCOV_EXCL_LINE */
1208 : }
1209 :
1210 : GEN
1211 270602775 : gaddsg(long x, GEN y)
1212 : {
1213 270602775 : long ty = typ(y);
1214 : GEN z;
1215 :
1216 270602775 : switch(ty)
1217 : {
1218 124843357 : case t_INT: return addsi(x,y);
1219 119204911 : case t_REAL: return addsr(x,y);
1220 12031 : case t_INTMOD:
1221 12031 : z = cgetg(3, t_INTMOD);
1222 12031 : return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1223 15247859 : case t_FRAC: z = cgetg(3,t_FRAC);
1224 15247859 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1225 15247859 : gel(z,2) = icopy(gel(y,2)); return z;
1226 8164474 : case t_COMPLEX:
1227 8164474 : z = cgetg(3, t_COMPLEX);
1228 8164474 : gel(z,1) = gaddsg(x, gel(y,1));
1229 8164474 : gel(z,2) = gcopy(gel(y,2)); return z;
1230 :
1231 3130143 : default: return gadd(stoi(x), y);
1232 : }
1233 : }
1234 :
1235 : GEN
1236 3168040 : gsubsg(long x, GEN y)
1237 : {
1238 : GEN z, a, b;
1239 : pari_sp av;
1240 :
1241 3168040 : switch(typ(y))
1242 : {
1243 276674 : case t_INT: return subsi(x,y);
1244 1225245 : case t_REAL: return subsr(x,y);
1245 56 : case t_INTMOD:
1246 56 : z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1247 56 : return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1248 732402 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1249 732402 : gel(z,1) = gc_INT((pari_sp)z, subii(mulis(b,x), a));
1250 732403 : gel(z,2) = icopy(gel(y,2)); return z;
1251 896549 : case t_COMPLEX:
1252 896549 : z = cgetg(3, t_COMPLEX);
1253 896549 : gel(z,1) = gsubsg(x, gel(y,1));
1254 896549 : gel(z,2) = gneg(gel(y,2)); return z;
1255 : }
1256 37114 : av = avma;
1257 37114 : return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
1258 : }
1259 :
1260 : /********************************************************************/
1261 : /** **/
1262 : /** SUBTRACTION **/
1263 : /** **/
1264 : /********************************************************************/
1265 :
1266 : GEN
1267 2848075709 : gsub(GEN x, GEN y)
1268 : {
1269 2848075709 : long tx = typ(x), ty = typ(y);
1270 : pari_sp av;
1271 : GEN z;
1272 2848075709 : if (tx == ty) switch(tx) /* shortcut to generic case */
1273 : {
1274 2068565872 : case t_INT: return subii(x,y);
1275 578280901 : case t_REAL: return subrr(x,y);
1276 1053735 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1277 1053735 : z = cgetg(3,t_INTMOD);
1278 1053735 : if (X==Y || equalii(X,Y))
1279 1053721 : return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1280 14 : gel(z,1) = gcdii(X,Y);
1281 14 : warn_coercion(X,Y,gel(z,1));
1282 14 : av = avma; p1 = subii(gel(x,2),gel(y,2));
1283 14 : gel(z,2) = gc_INT(av, modii(p1, gel(z,1))); return z;
1284 : }
1285 5071148 : case t_FRAC: return addsub_frac(x,y, subii);
1286 102062433 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1287 102098337 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1288 102007979 : if (isintzero(gel(z,2)))
1289 : {
1290 21301 : set_avma((pari_sp)(z+3));
1291 21301 : return gsub(gel(x,1),gel(y,1));
1292 : }
1293 101989764 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1294 101974718 : return z;
1295 5438381 : case t_PADIC:
1296 5438381 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1297 5438381 : return addsub_pp(x,y, subii);
1298 168 : case t_QUAD: z = cgetg(4,t_QUAD);
1299 168 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1300 168 : gel(z,1) = ZX_copy(gel(x,1));
1301 168 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1302 168 : gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1303 850444 : case t_POLMOD:
1304 850444 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1305 850409 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1306 35 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1307 240376 : case t_FFELT: return FF_sub(x,y);
1308 18419209 : case t_POL: {
1309 18419209 : long vx = varn(x);
1310 18419209 : long vy = varn(y);
1311 18419209 : if (vx != vy) {
1312 1146754 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1313 84813 : else return Rg_RgX_sub(x, y);
1314 : }
1315 17272455 : return RgX_sub(x, y);
1316 : }
1317 299567 : case t_VEC:
1318 299567 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1319 299567 : return RgV_sub(x,y);
1320 3568664 : case t_COL:
1321 3568664 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1322 3568664 : return RgC_sub(x,y);
1323 252131 : case t_MAT: {
1324 252131 : long lx = lg(x);
1325 252131 : if (lg(y) != lx) pari_err_OP("+",x,y);
1326 252133 : if (lx == 1) return cgetg(1, t_MAT);
1327 172505 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1328 172505 : return RgM_sub(x,y);
1329 : }
1330 2448296 : case t_RFRAC: case t_SER: break;
1331 :
1332 0 : default: pari_err_TYPE2("+",x,y);
1333 : }
1334 65609357 : av = avma;
1335 65609357 : return gerepileupto(av, gadd(x,gneg_i(y)));
1336 : }
1337 :
1338 : /********************************************************************/
1339 : /** **/
1340 : /** MULTIPLICATION **/
1341 : /** **/
1342 : /********************************************************************/
1343 : static GEN
1344 316204 : mul_ser_scal(GEN x, GEN t)
1345 : {
1346 316204 : if (isexactzero(t)) return gmul(Rg_get_0(x), t);
1347 312956 : if (isint1(t)) return gcopy(x);
1348 262388 : if (ser_isexactzero(x))
1349 : {
1350 378 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gmul(gel(x,2), t), varn(x), 1);
1351 378 : setvalser(z, valser(x)); return z;
1352 : }
1353 3442418 : pari_APPLY_ser(gmul(gel(x,i), t));
1354 : }
1355 : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1356 : * [n/d a valid RFRAC] */
1357 : static GEN
1358 10416663 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1359 : {
1360 10416663 : pari_sp av = avma;
1361 : GEN z;
1362 :
1363 10416663 : switch(typ(x))
1364 : {
1365 21 : case t_PADIC:
1366 21 : n = gmul(n, x);
1367 21 : d = gcvtop(d, padic_p(x), signe(padic_u(x))? precp(x): 1);
1368 21 : return gerepileupto(av, gdiv(n,d));
1369 :
1370 623 : case t_INTMOD: case t_POLMOD:
1371 623 : n = gmul(n, x);
1372 623 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1373 623 : return gerepileupto(av, gdiv(n,d));
1374 : }
1375 10416019 : z = gred_rfrac2(x, d);
1376 10416019 : n = simplify_shallow(n);
1377 10416019 : if (typ(z) == t_RFRAC)
1378 : {
1379 7897157 : n = gmul(gel(z,1), n);
1380 7897157 : d = gel(z,2);
1381 7897157 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1382 0 : z = RgX_Rg_div(n, d);
1383 : else
1384 7897157 : z = gred_rfrac_simple(n, d);
1385 : }
1386 : else
1387 2518862 : z = gmul(z, n);
1388 10416019 : return gerepileupto(av, z);
1389 : }
1390 : static GEN
1391 150273552 : mul_scal(GEN y, GEN x, long ty)
1392 : {
1393 150273552 : switch(ty)
1394 : {
1395 141379150 : case t_POL:
1396 141379150 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1397 129460232 : return RgX_Rg_mul(y, x);
1398 308063 : case t_SER: return mul_ser_scal(y, x);
1399 8586326 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1400 14 : case t_QFB:
1401 14 : if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1402 : }
1403 12 : pari_err_TYPE2("*",x,y);
1404 : return NULL; /* LCOV_EXCL_LINE */
1405 : }
1406 : static GEN
1407 7461129 : mul_self_scal(GEN x, GEN y)
1408 633925443 : { pari_APPLY_same(gmul(y,gel(x,i))); }
1409 :
1410 : static GEN
1411 160872 : mul_gen_rfrac(GEN X, GEN Y)
1412 : {
1413 160872 : GEN y1 = gel(Y,1), y2 = gel(Y,2);
1414 160872 : long vx = gvar(X), vy = varn(y2);
1415 166633 : return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1416 5761 : gred_rfrac_simple(gmul(y1,X), y2);
1417 : }
1418 : /* (x1/x2) * (y1/y2) */
1419 : static GEN
1420 7907758 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1421 : {
1422 : GEN z, X, Y;
1423 7907758 : pari_sp av = avma;
1424 :
1425 7907758 : X = gred_rfrac2(x1, y2);
1426 7907758 : Y = gred_rfrac2(y1, x2);
1427 7907758 : if (typ(X) == t_RFRAC)
1428 : {
1429 6628606 : if (typ(Y) == t_RFRAC) {
1430 6563009 : x1 = gel(X,1);
1431 6563009 : x2 = gel(X,2);
1432 6563009 : y1 = gel(Y,1);
1433 6563009 : y2 = gel(Y,2);
1434 6563009 : z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1435 : } else
1436 65597 : z = mul_gen_rfrac(Y, X);
1437 : }
1438 1279152 : else if (typ(Y) == t_RFRAC)
1439 95275 : z = mul_gen_rfrac(X, Y);
1440 : else
1441 1183877 : z = gmul(X, Y);
1442 7907758 : return gerepileupto(av, z);
1443 : }
1444 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1445 : static GEN
1446 262512 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1447 : {
1448 262512 : pari_sp av = avma;
1449 262512 : GEN X = gred_rfrac2(x1, y2);
1450 262512 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1451 : {
1452 255932 : x2 = RgX_mul(gel(X,2), x2);
1453 255932 : x1 = gel(X,1);
1454 : }
1455 : else
1456 6580 : x1 = X;
1457 262512 : return gerepileupto(av, gred_rfrac_simple(x1, x2));
1458 : }
1459 :
1460 : /* Mod(y, Y) * x, assuming x scalar */
1461 : static GEN
1462 3428782 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1463 3428782 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1464 :
1465 : /* cf mulqq */
1466 : static GEN
1467 5860762 : quad_polmod_mul(GEN P, GEN x, GEN y)
1468 : {
1469 5860762 : GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
1470 5860762 : pari_sp tetpil, av = avma;
1471 5860762 : T[1] = x[1];
1472 5860762 : p2 = gmul(gel(x,2), gel(y,2));
1473 5860762 : p3 = gmul(gel(x,3), gel(y,3));
1474 5860762 : p1 = gmul(gneg_i(c),p3);
1475 : /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
1476 5860762 : if (typ(b) == t_INT)
1477 : {
1478 5860741 : if (signe(b))
1479 : {
1480 4284989 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1481 4284989 : if (is_pm1(b))
1482 : {
1483 4284331 : if (signe(b) > 0) p3 = gneg(p3);
1484 : }
1485 : else
1486 658 : p3 = gmul(negi(b), p3);
1487 : }
1488 : else
1489 : {
1490 1575752 : p3 = gmul(gel(x,2),gel(y,3));
1491 1575752 : p4 = gmul(gel(x,3),gel(y,2));
1492 : }
1493 : }
1494 : else
1495 : {
1496 21 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1497 21 : p3 = gmul(gneg_i(b), p3);
1498 : }
1499 5860762 : tetpil = avma;
1500 5860762 : gel(T,2) = gadd(p2, p1);
1501 5860762 : gel(T,3) = gadd(p4, p3);
1502 5860762 : gerepilecoeffssp(av,tetpil,T+2,2);
1503 5860762 : return normalizepol_lg(T,4);
1504 : }
1505 : /* Mod(x,T) * Mod(y,T) */
1506 : static GEN
1507 8565249 : mul_polmod_same(GEN T, GEN x, GEN y)
1508 : {
1509 8565249 : GEN z = cgetg(3,t_POLMOD), a;
1510 8565249 : long v = varn(T), lx = lg(x), ly = lg(y);
1511 8565249 : gel(z,1) = RgX_copy(T);
1512 : /* x * y mod T optimised */
1513 8565249 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1514 7894305 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1515 1588179 : a = gmul(x, y);
1516 : else
1517 : {
1518 6977070 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1519 5855484 : a = quad_polmod_mul(T, x, y);
1520 : else
1521 1121586 : a = RgXQ_mul(x, y, gel(z,1));
1522 : }
1523 8565249 : gel(z,2) = a; return z;
1524 : }
1525 : static GEN
1526 484275 : sqr_polmod(GEN T, GEN x)
1527 : {
1528 484275 : GEN a, z = cgetg(3,t_POLMOD);
1529 484275 : gel(z,1) = RgX_copy(T);
1530 484275 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1531 93885 : a = gsqr(x);
1532 : else
1533 : {
1534 390390 : pari_sp av = avma;
1535 390390 : a = RgXQ_sqr(x, gel(z,1));
1536 390390 : a = gerepileupto(av, a);
1537 : }
1538 484275 : gel(z,2) = a; return z;
1539 : }
1540 : /* Mod(x,X) * Mod(y,Y) */
1541 : static GEN
1542 2674917 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1543 : {
1544 2674917 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1545 2674917 : long vx = varn(X), vy = varn(Y);
1546 2674917 : GEN z = cgetg(3,t_POLMOD);
1547 :
1548 2674917 : if (vx==vy) {
1549 : pari_sp av;
1550 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
1551 14 : warn_coercion(X,Y,gel(z,1));
1552 14 : gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
1553 14 : return z;
1554 : }
1555 2674903 : if (varncmp(vx, vy) < 0)
1556 410662 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1557 : else
1558 2264241 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1559 2674903 : gel(z,2) = gmul(x, y); return z;
1560 : }
1561 :
1562 : #if 0 /* used by 3M only */
1563 : /* set z = x+y and return 1 if x,y have the same sign
1564 : * set z = x-y and return 0 otherwise */
1565 : static int
1566 : did_add(GEN x, GEN y, GEN *z)
1567 : {
1568 : long tx = typ(x), ty = typ(y);
1569 : if (tx == ty) switch(tx)
1570 : {
1571 : case t_INT: *z = addii(x,y); return 1;
1572 : case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1573 : case t_REAL:
1574 : if (signe(x) == -signe(y))
1575 : { *z = subrr(x,y); return 0; }
1576 : else
1577 : { *z = addrr(x,y); return 1; }
1578 : }
1579 : if (tx == t_REAL) switch(ty)
1580 : {
1581 : case t_INT:
1582 : if (signe(x) == -signe(y))
1583 : { *z = subri(x,y); return 0; }
1584 : else
1585 : { *z = addri(x,y); return 1; }
1586 : case t_FRAC:
1587 : if (signe(x) == -signe(gel(y,1)))
1588 : { *z = gsub(x,y); return 0; }
1589 : else
1590 : { *z = gadd(x,y); return 1; }
1591 : }
1592 : else if (ty == t_REAL) switch(tx)
1593 : {
1594 : case t_INT:
1595 : if (signe(x) == -signe(y))
1596 : { *z = subir(x,y); return 0; }
1597 : else
1598 : { *z = addir(x,y); return 1; }
1599 : case t_FRAC:
1600 : if (signe(gel(x,1)) == -signe(y))
1601 : { *z = gsub(x,y); return 0; }
1602 : else
1603 : { *z = gadd(x,y); return 1; }
1604 : }
1605 : *z = gadd(x,y); return 1;
1606 : }
1607 : #endif
1608 : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1609 : static GEN
1610 11664518 : mulcIR(GEN x, GEN y)
1611 : {
1612 11664518 : GEN z = cgetg(3,t_COMPLEX);
1613 11664671 : pari_sp av = avma;
1614 11664671 : gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
1615 11664939 : gel(z,2) = gmul(y, gel(x,1));
1616 11664800 : return z;
1617 :
1618 : }
1619 : /* x,y COMPLEX */
1620 : static GEN
1621 273943998 : mulcc(GEN x, GEN y)
1622 : {
1623 273943998 : GEN xr = gel(x,1), xi = gel(x,2);
1624 273943998 : GEN yr = gel(y,1), yi = gel(y,2);
1625 : GEN p1, p2, p3, p4, z;
1626 : pari_sp tetpil, av;
1627 :
1628 273943998 : if (isintzero(xr))
1629 : {
1630 15646141 : if (isintzero(yr)) {
1631 7128913 : av = avma;
1632 7128913 : return gerepileupto(av, gneg(gmul(xi,yi)));
1633 : }
1634 8517046 : return mulcIR(y, xi);
1635 : }
1636 258286077 : if (isintzero(yr)) return mulcIR(x, yi);
1637 :
1638 255127866 : z = cgetg(3,t_COMPLEX); av = avma;
1639 : #if 0
1640 : /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1641 : * e.g. xr + xi if exponents differ */
1642 : if (did_add(xr, xi, &p3))
1643 : {
1644 : if (did_add(yr, yi, &p4)) {
1645 : /* R = xr*yr - xi*yi
1646 : * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1647 : p1 = gmul(xr,yr);
1648 : p2 = gmul(xi,yi); p2 = gneg(p2);
1649 : p3 = gmul(p3, p4);
1650 : p4 = gsub(p2, p1);
1651 : } else {
1652 : /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1653 : * I = xr*yi + xi*yr */
1654 : p1 = gmul(p3,p4);
1655 : p3 = gmul(xr,yi);
1656 : p4 = gmul(xi,yr);
1657 : p2 = gsub(p3, p4);
1658 : }
1659 : } else {
1660 : if (did_add(yr, yi, &p4)) {
1661 : /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1662 : * I = xr*yi +xi*yr */
1663 : p1 = gmul(p3,p4);
1664 : p3 = gmul(xr,yi);
1665 : p4 = gmul(xi,yr);
1666 : p2 = gsub(p4, p3);
1667 : } else {
1668 : /* R = xr*yr - xi*yi
1669 : * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1670 : p3 = gneg( gmul(p3, p4) );
1671 : p1 = gmul(xr,yr);
1672 : p2 = gmul(xi,yi);
1673 : p4 = gadd(p1, p2);
1674 :
1675 : p2 = gneg(p2);
1676 : }
1677 : }
1678 : tetpil = avma;
1679 : gel(z,1) = gadd(p1,p2);
1680 : gel(z,2) = gadd(p3,p4);
1681 : #else
1682 255130432 : if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1683 : { /* 3M formula */
1684 559356 : p3 = addii(xr,xi);
1685 559356 : p4 = addii(yr,yi);
1686 559356 : p1 = mulii(xr,yr);
1687 559356 : p2 = mulii(xi,yi);
1688 559356 : p3 = mulii(p3,p4);
1689 559356 : p4 = addii(p2,p1);
1690 559356 : tetpil = avma;
1691 559356 : gel(z,1) = subii(p1,p2);
1692 559356 : gel(z,2) = subii(p3,p4);
1693 559356 : if (!signe(gel(z,2)))
1694 113260 : return gc_INT((pari_sp)(z+3), gel(z,1));
1695 : }
1696 : else
1697 : { /* naive 4M formula: avoid all loss of accuracy */
1698 254571076 : p1 = gmul(xr,yr);
1699 254528166 : p2 = gmul(xi,yi);
1700 254529636 : p3 = gmul(xr,yi);
1701 254532437 : p4 = gmul(xi,yr);
1702 254541361 : tetpil = avma;
1703 254541361 : gel(z,1) = gsub(p1,p2);
1704 254395749 : gel(z,2) = gadd(p3,p4);
1705 254409079 : if (isintzero(gel(z,2)))
1706 : {
1707 50575 : cgiv(gel(z,2));
1708 50575 : return gerepileupto((pari_sp)(z+3), gel(z,1));
1709 : }
1710 : }
1711 : #endif
1712 :
1713 254802490 : gerepilecoeffssp(av,tetpil, z+1,2); return z;
1714 : }
1715 : /* x,y PADIC */
1716 : static GEN
1717 17503103 : mulpp(GEN x, GEN y) {
1718 17503103 : GEN pd, p = padic_p(x), ux = padic_u(x), uy = padic_u(y);
1719 17503103 : long e = valp(x) + valp(y), dx, dy;
1720 :
1721 17503103 : if (!equalii(p, padic_p(y))) pari_err_OP("*",x,y);
1722 17503056 : if (!signe(ux) || !signe(uy)) return zeropadic(p, e);
1723 16930305 : dx = precp(x); dy = precp(y);
1724 16930305 : if (dx > dy) { pd = padic_pd(y); dx = dy; } else pd = padic_pd(x);
1725 16930305 : retmkpadic(Fp_mul(ux, uy, pd), icopy(p), icopy(pd), e, dx);
1726 : }
1727 : /* x,y QUAD */
1728 : static GEN
1729 1106 : mulqq(GEN x, GEN y) {
1730 1106 : GEN z = cgetg(4,t_QUAD);
1731 1106 : GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
1732 : pari_sp av, tetpil;
1733 1106 : if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
1734 :
1735 1106 : gel(z,1) = ZX_copy(P); av = avma;
1736 1106 : p2 = gmul(gel(x,2),gel(y,2));
1737 1106 : p3 = gmul(gel(x,3),gel(y,3));
1738 1106 : p1 = gmul(gneg_i(c),p3);
1739 :
1740 1106 : if (signe(b))
1741 987 : p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
1742 : else
1743 : {
1744 119 : p3 = gmul(gel(x,2),gel(y,3));
1745 119 : p4 = gmul(gel(x,3),gel(y,2));
1746 : }
1747 1106 : tetpil = avma;
1748 1106 : gel(z,2) = gadd(p2,p1);
1749 1106 : gel(z,3) = gadd(p4,p3);
1750 1106 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
1751 : }
1752 :
1753 : GEN
1754 15648965 : mulcxI(GEN x)
1755 : {
1756 : GEN z;
1757 15648965 : switch(typ(x))
1758 : {
1759 1006073 : case t_INT:
1760 1006073 : if (!signe(x)) return gen_0;
1761 : case t_REAL: case t_FRAC:
1762 2631948 : return mkcomplex(gen_0, x);
1763 13016786 : case t_COMPLEX:
1764 13016786 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1765 13003313 : z = cgetg(3,t_COMPLEX);
1766 13003635 : gel(z,1) = gneg(gel(x,2));
1767 13004268 : gel(z,2) = gel(x,1); return z;
1768 56 : default:
1769 56 : return gmul(gen_I(), x);
1770 : }
1771 : }
1772 : GEN
1773 3079001 : mulcxmI(GEN x)
1774 : {
1775 : GEN z;
1776 3079001 : switch(typ(x))
1777 : {
1778 1106 : case t_INT:
1779 1106 : if (!signe(x)) return gen_0;
1780 : case t_REAL: case t_FRAC:
1781 117193 : return mkcomplex(gen_0, gneg(x));
1782 2758904 : case t_COMPLEX:
1783 2758904 : if (isintzero(gel(x,1))) return gel(x,2);
1784 2720706 : z = cgetg(3,t_COMPLEX);
1785 2720679 : gel(z,1) = gel(x,2);
1786 2720679 : gel(z,2) = gneg(gel(x,1)); return z;
1787 202904 : default:
1788 202904 : return gmul(mkcomplex(gen_0, gen_m1), x);
1789 : }
1790 : }
1791 : /* x * I^k */
1792 : GEN
1793 5523587 : mulcxpowIs(GEN x, long k)
1794 : {
1795 5523587 : switch (k & 3)
1796 : {
1797 1369148 : case 1: return mulcxI(x);
1798 1362180 : case 2: return gneg(x);
1799 1363811 : case 3: return mulcxmI(x);
1800 : }
1801 1428448 : return x;
1802 : }
1803 :
1804 : /* caller will assume l > 2 later */
1805 : static GEN
1806 5706235 : init_ser(long l, long v, long e)
1807 : {
1808 5706235 : GEN z = cgetg(l, t_SER);
1809 5706235 : z[1] = evalvalser(e) | evalvarn(v) | evalsigne(1); return z;
1810 : }
1811 :
1812 : /* fill in coefficients of t_SER z from coeffs of t_POL y */
1813 : static GEN
1814 5693888 : fill_ser(GEN z, GEN y)
1815 : {
1816 5693888 : long i, lx = lg(z), ly = lg(y); /* lx > 2 */
1817 5693888 : if (ly >= lx) {
1818 25570501 : for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
1819 : } else {
1820 335415 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1821 239573 : for ( ; i < lx; i++) gel(z,i) = gen_0;
1822 : }
1823 : /* dangerous case (already normalized), don't use normalizeser */
1824 5693888 : if (ly == 3 && !signe(y)) { setsigne(z, 0); return z; }
1825 5693790 : return normalizeser(z);
1826 : }
1827 : /* assume typ(x) = t_VEC */
1828 : static int
1829 98 : is_ext_qfr(GEN x)
1830 91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
1831 189 : && typ(gel(x,2)) == t_REAL; }
1832 :
1833 : GEN
1834 8754531763 : gmul(GEN x, GEN y)
1835 : {
1836 : long tx, ty, lx, ly, vx, vy, i, l;
1837 : pari_sp av, tetpil;
1838 : GEN z, p1;
1839 :
1840 8754531763 : if (x == y) return gsqr(x);
1841 7820779015 : tx = typ(x); ty = typ(y);
1842 7820779015 : if (tx == ty) switch(tx)
1843 : {
1844 3663705431 : case t_INT: return mulii(x,y);
1845 1869600224 : case t_REAL: return mulrr(x,y);
1846 1545015 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1847 1545015 : z = cgetg(3,t_INTMOD);
1848 1545015 : if (X==Y || equalii(X,Y))
1849 1544980 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1850 35 : gel(z,1) = gcdii(X,Y);
1851 35 : warn_coercion(X,Y,gel(z,1));
1852 35 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1853 35 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
1854 : }
1855 23261392 : case t_FRAC:
1856 : {
1857 23261392 : GEN x1 = gel(x,1), x2 = gel(x,2);
1858 23261392 : GEN y1 = gel(y,1), y2 = gel(y,2);
1859 23261392 : z=cgetg(3,t_FRAC);
1860 23261392 : p1 = gcdii(x1, y2);
1861 23261388 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1862 23261388 : p1 = gcdii(x2, y1);
1863 23261390 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1864 23261390 : tetpil = avma;
1865 23261390 : gel(z,2) = mulii(x2,y2);
1866 23261389 : gel(z,1) = mulii(x1,y1);
1867 23261392 : fix_frac_if_int_GC(z,tetpil); return z;
1868 : }
1869 264941791 : case t_COMPLEX: return mulcc(x, y);
1870 12511724 : case t_PADIC: return mulpp(x, y);
1871 882 : case t_QUAD: return mulqq(x, y);
1872 13693856 : case t_FFELT: return FF_mul(x, y);
1873 11004294 : case t_POLMOD:
1874 11004294 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1875 8329377 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1876 2674917 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1877 128360480 : case t_POL:
1878 128360480 : vx = varn(x);
1879 128360480 : vy = varn(y);
1880 128360480 : if (vx != vy) {
1881 45999793 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1882 26874702 : else return RgX_Rg_mul(y, x);
1883 : }
1884 82360687 : return RgX_mul(x, y);
1885 :
1886 4411812 : case t_SER: {
1887 4411812 : vx = varn(x);
1888 4411812 : vy = varn(y);
1889 4411812 : if (vx != vy) {
1890 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1891 3675 : return mul_ser_scal(y, x);
1892 : }
1893 4408137 : lx = minss(lg(x), lg(y));
1894 4408137 : if (lx == 2) return zeroser(vx, valser(x)+valser(y));
1895 4408123 : av = avma; z = init_ser(lx, vx, valser(x)+valser(y));
1896 4408123 : x = ser2pol_i(x, lx);
1897 4408123 : y = ser2pol_i(y, lx);
1898 4408123 : y = RgXn_mul(x, y, lx-2);
1899 4408123 : return gc_GEN(av, fill_ser(z,y));
1900 : }
1901 42 : case t_VEC:
1902 42 : if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
1903 : /* fall through, handle extended t_QFB */
1904 45632 : case t_QFB: return qfbcomp(x,y);
1905 6721025 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1906 4114815 : case t_MAT: return RgM_mul(x, y);
1907 :
1908 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1909 1631 : z = cgetg_copy(x, &l);
1910 1631 : if (l != lg(y)) break;
1911 17325 : for (i=1; i<l; i++)
1912 : {
1913 15694 : long yi = y[i];
1914 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1915 15694 : z[i] = x[yi];
1916 : }
1917 1631 : return z;
1918 :
1919 0 : default:
1920 0 : pari_err_TYPE2("*",x,y);
1921 : }
1922 : /* tx != ty */
1923 1819386562 : if (is_const_t(ty) && is_const_t(tx)) {
1924 1660616318 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1925 1660616318 : switch(tx) {
1926 1531694841 : case t_INT:
1927 : switch(ty)
1928 : {
1929 987846169 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1930 1290437 : case t_INTMOD:
1931 1290437 : z = cgetg(3, t_INTMOD);
1932 1290436 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1933 105930077 : case t_FRAC:
1934 105930077 : if (!signe(x)) return gen_0;
1935 88882841 : z=cgetg(3,t_FRAC);
1936 88883087 : p1 = gcdii(x,gel(y,2));
1937 88882161 : if (equali1(p1))
1938 : {
1939 49186913 : set_avma((pari_sp)z);
1940 49186899 : gel(z,2) = icopy(gel(y,2));
1941 49186991 : gel(z,1) = mulii(gel(y,1), x);
1942 : }
1943 : else
1944 : {
1945 39695579 : x = diviiexact(x,p1); tetpil = avma;
1946 39695046 : gel(z,2) = diviiexact(gel(y,2), p1);
1947 39694730 : gel(z,1) = mulii(gel(y,1), x);
1948 39695189 : fix_frac_if_int_GC(z,tetpil);
1949 : }
1950 88882730 : return z;
1951 433756256 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1952 4196889 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1953 1659 : case t_QUAD: return mulRq(x,y);
1954 1607192 : case t_FFELT: return FF_Z_mul(y,x);
1955 : }
1956 :
1957 : case t_REAL:
1958 121254751 : switch(ty)
1959 : {
1960 37585302 : case t_FRAC: return mulrfrac(x, y);
1961 83669407 : case t_COMPLEX: return mulRc(x, y);
1962 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1963 21 : default: pari_err_TYPE2("*",x,y);
1964 : }
1965 :
1966 8036 : case t_INTMOD:
1967 : switch(ty)
1968 : {
1969 7133 : case t_FRAC: { GEN X = gel(x,1);
1970 7133 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1971 7133 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1972 : }
1973 49 : case t_COMPLEX: return mulRc_direct(x,y);
1974 259 : case t_PADIC: { GEN X = gel(x,1);
1975 259 : z = cgetg(3, t_INTMOD);
1976 259 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1977 : }
1978 63 : case t_QUAD: return mulRq(x, y);
1979 532 : case t_FFELT:
1980 532 : if (!equalii(gel(x,1),FF_p_i(y)))
1981 0 : pari_err_OP("*",x,y);
1982 532 : return FF_Z_mul(y,gel(x,2));
1983 : }
1984 :
1985 : case t_FRAC:
1986 : switch(ty)
1987 : {
1988 5272610 : case t_COMPLEX: return mulRc(x, y);
1989 2582756 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
1990 343 : case t_QUAD: return mulRq(x, y);
1991 2051 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
1992 : }
1993 :
1994 : case t_FFELT:
1995 35 : pari_err_TYPE2("*",x,y);
1996 :
1997 21 : case t_COMPLEX:
1998 : switch(ty)
1999 : {
2000 14 : case t_PADIC:
2001 14 : return Zp_nosquare_m1(padic_p(y))? mulRc(y, x): mulTp(x, y);
2002 7 : case t_QUAD:
2003 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
2004 7 : return mulqf(y, x, lx);
2005 : }
2006 :
2007 : case t_PADIC: /* ty == t_QUAD */
2008 28 : return (kro_quad(y,padic_p(x))== -1)? mulRq(x, y): mulTp(y, x);
2009 : }
2010 : }
2011 :
2012 160752851 : if (is_matvec_t(ty)) switch(tx)
2013 : {
2014 122530 : case t_VEC:
2015 : switch(ty) {
2016 15540 : case t_COL: return RgV_RgC_mul(x,y);
2017 106990 : case t_MAT: return RgV_RgM_mul(x,y);
2018 : }
2019 0 : break;
2020 1687 : case t_COL:
2021 : switch(ty) {
2022 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2023 0 : case t_MAT: return RgC_RgM_mul(x,y);
2024 : }
2025 0 : break;
2026 819645 : case t_MAT:
2027 : switch(ty) {
2028 0 : case t_VEC: return RgM_RgV_mul(x,y);
2029 819646 : case t_COL: return RgM_RgC_mul(x,y);
2030 : }
2031 : default:
2032 5403253 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
2033 5403258 : return mul_self_scal(y, x);
2034 : }
2035 157675365 : if (is_matvec_t(tx))
2036 : {
2037 2057580 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2038 2057825 : return mul_self_scal(x, y);
2039 : }
2040 155618324 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2041 : /* tx < ty, !ismatvec(x and y) */
2042 :
2043 155618324 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2044 3405689 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2045 152212635 : if (is_scalar_t(tx)) {
2046 148957238 : if (tx == t_POLMOD) {
2047 3110312 : vx = varn(gel(x,1));
2048 3110312 : vy = gvar(y);
2049 3110312 : if (vx != vy) {
2050 2875000 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2051 23093 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2052 : }
2053 : /* error if ty == t_SER */
2054 235312 : av = avma; y = gmod(y, gel(x,1));
2055 235305 : return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2056 : }
2057 145846926 : return mul_scal(y, x, ty);
2058 : }
2059 :
2060 : /* x and y are not scalars, nor matvec */
2061 3255292 : vx = gvar(x);
2062 3255327 : vy = gvar(y);
2063 3255327 : if (vx != vy) /* x or y is treated as a scalar */
2064 2781167 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2065 2781167 : : mul_scal(y, x, ty);
2066 : /* vx = vy */
2067 1835949 : switch(tx)
2068 : {
2069 1835914 : case t_POL:
2070 : switch (ty)
2071 : {
2072 6741 : case t_SER:
2073 : {
2074 : long v;
2075 6741 : av = avma; v = RgX_valrem(x, &x);
2076 6741 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2077 6727 : v += valser(y); ly = lg(y);
2078 6727 : if (ly == 2) { set_avma(av); return zeroser(vx, v); }
2079 6496 : if (degpol(x))
2080 : {
2081 2030 : z = init_ser(ly, vx, v);
2082 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2083 2030 : return gc_GEN(av, fill_ser(z, x));
2084 : }
2085 : /* take advantage of x = c*t^v */
2086 4466 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2087 4466 : setvalser(y, v); return y;
2088 : }
2089 :
2090 1829173 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2091 : }
2092 0 : break;
2093 :
2094 35 : case t_SER:
2095 : switch (ty)
2096 : {
2097 35 : case t_RFRAC:
2098 35 : av = avma;
2099 35 : return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2100 : }
2101 0 : break;
2102 : }
2103 0 : pari_err_TYPE2("*",x,y);
2104 : return NULL; /* LCOV_EXCL_LINE */
2105 : }
2106 :
2107 : /* return a nonnormalized result */
2108 : GEN
2109 112731 : sqr_ser_part(GEN x, long l1, long l2)
2110 : {
2111 : long i, j, l;
2112 : pari_sp av;
2113 : GEN Z, z, p1, p2;
2114 : long mi;
2115 112731 : if (l2 < l1) return zeroser(varn(x), 2*valser(x));
2116 112717 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2117 112717 : Z = cgetg(l2-l1+3, t_SER);
2118 112717 : Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
2119 112717 : z = Z + 2-l1;
2120 112717 : x += 2; mi = 0;
2121 425390 : for (i=0; i<l1; i++)
2122 : {
2123 312673 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2124 : }
2125 :
2126 720212 : for (i=l1; i<=l2; i++)
2127 : {
2128 607495 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2129 607495 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2130 1256134 : for (j=i-mi; j<=minss(l,mi); j++)
2131 648639 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2132 607495 : p1 = gshift(p1,1);
2133 607495 : if ((i&1) == 0 && p2[i>>1])
2134 93585 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2135 607495 : gel(z,i) = gerepileupto(av,p1);
2136 : }
2137 112717 : return Z;
2138 : }
2139 :
2140 : GEN
2141 1175201223 : gsqr(GEN x)
2142 : {
2143 : long i, lx;
2144 : pari_sp av, tetpil;
2145 : GEN z, p1, p2, p3, p4;
2146 :
2147 1175201223 : switch(typ(x))
2148 : {
2149 962560896 : case t_INT: return sqri(x);
2150 195333370 : case t_REAL: return sqrr(x);
2151 142489 : case t_INTMOD: { GEN X = gel(x,1);
2152 142489 : z = cgetg(3,t_INTMOD);
2153 142489 : gel(z,2) = gc_INT((pari_sp)z, remii(sqri(gel(x,2)), X));
2154 142489 : gel(z,1) = icopy(X); return z;
2155 : }
2156 4136070 : case t_FRAC: return sqrfrac(x);
2157 :
2158 8117704 : case t_COMPLEX:
2159 8117704 : if (isintzero(gel(x,1))) {
2160 238697 : av = avma;
2161 238697 : return gerepileupto(av, gneg(gsqr(gel(x,2))));
2162 : }
2163 7879006 : z = cgetg(3,t_COMPLEX); av = avma;
2164 7878995 : p1 = gadd(gel(x,1),gel(x,2));
2165 7878811 : p2 = gsub(gel(x,1), gel(x,2));
2166 7878803 : p3 = gmul(gel(x,1),gel(x,2));
2167 7878933 : tetpil = avma;
2168 7878933 : gel(z,1) = gmul(p1,p2);
2169 7878994 : gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
2170 :
2171 4795 : case t_PADIC:
2172 : {
2173 4795 : GEN u = padic_u(x), p = padic_p(x), pd = padic_pd(x);
2174 4795 : long v2 = 2*valp(x), d = precp(x);
2175 4795 : if (!signe(u)) { x = gcopy(x); setvalp(x, v2); return x; }
2176 4739 : if (!absequaliu(p, 2))
2177 3206 : retmkpadic(Fp_sqr(u, pd), icopy(p), icopy(pd), v2, d);
2178 : /* p = 2*/
2179 1533 : if (d == 1) /* (1 + O(2))^2 = 1 + O(2^3) */
2180 7 : retmkpadic(gen_1, gen_2, utoipos(8), v2, 3);
2181 1526 : retmkpadic(remi2n(sqri(u), d + 1), gen_2, int2n(d + 1), v2, d + 1);
2182 : }
2183 70 : case t_QUAD: z = cgetg(4,t_QUAD);
2184 70 : p1 = gel(x,1);
2185 70 : gel(z,1) = ZX_copy(p1); av = avma;
2186 70 : p2 = gsqr(gel(x,2));
2187 70 : p3 = gsqr(gel(x,3));
2188 70 : p4 = gmul(gneg_i(gel(p1,2)),p3);
2189 :
2190 70 : if (gequal0(gel(p1,3)))
2191 : {
2192 7 : tetpil = avma;
2193 7 : gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
2194 7 : av = avma;
2195 7 : p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
2196 7 : gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
2197 : }
2198 :
2199 63 : p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
2200 63 : tetpil = avma;
2201 63 : gel(z,2) = gadd(p2,p4);
2202 63 : gel(z,3) = gadd(p1,p3);
2203 63 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
2204 :
2205 484275 : case t_POLMOD:
2206 484275 : return sqr_polmod(gel(x,1), gel(x,2));
2207 :
2208 2967120 : case t_FFELT: return FF_sqr(x);
2209 :
2210 1352305 : case t_POL: return RgX_sqr(x);
2211 :
2212 35014 : case t_SER:
2213 35014 : lx = lg(x);
2214 35014 : if (ser_isexactzero(x)) {
2215 21 : GEN z = gcopy(x);
2216 21 : setvalser(z, 2*valser(x));
2217 21 : return z;
2218 : }
2219 34993 : if (lx < 40)
2220 34734 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2221 : else
2222 : {
2223 259 : pari_sp av = avma;
2224 259 : GEN z = init_ser(lx, varn(x), 2*valser(x));
2225 259 : x = ser2pol_i(x, lx);
2226 259 : x = RgXn_sqr(x, lx-2);
2227 259 : return gc_GEN(av, fill_ser(z,x));
2228 : }
2229 :
2230 63 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2231 63 : gel(z,1) = gsqr(gel(x,1));
2232 63 : gel(z,2) = gsqr(gel(x,2)); return z;
2233 :
2234 1169 : case t_MAT: return RgM_sqr(x);
2235 28 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
2236 : /* fall through handle extended t_QFB */
2237 80245 : case t_QFB: return qfbsqr(x);
2238 658 : case t_VECSMALL:
2239 658 : z = cgetg_copy(x, &lx);
2240 16289 : for (i=1; i<lx; i++)
2241 : {
2242 15631 : long xi = x[i];
2243 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2244 15631 : z[i] = x[xi];
2245 : }
2246 658 : return z;
2247 : }
2248 7 : pari_err_TYPE2("*",x,x);
2249 : return NULL; /* LCOV_EXCL_LINE */
2250 : }
2251 :
2252 : /********************************************************************/
2253 : /** **/
2254 : /** DIVISION **/
2255 : /** **/
2256 : /********************************************************************/
2257 : static GEN
2258 13162 : div_rfrac_scal(GEN x, GEN y)
2259 : {
2260 13162 : pari_sp av = avma;
2261 13162 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2262 13162 : return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
2263 : }
2264 : static GEN
2265 37515 : div_scal_rfrac(GEN x, GEN y)
2266 : {
2267 37515 : GEN y1 = gel(y,1), y2 = gel(y,2);
2268 37515 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2269 : {
2270 14 : if (degpol(y1))
2271 : {
2272 14 : pari_sp av = avma;
2273 14 : GEN _1 = Rg_get_1(x);
2274 14 : if (transtype(_1)) y1 = gmul(y1, _1);
2275 14 : return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
2276 : }
2277 0 : y1 = gel(y1,2);
2278 : }
2279 37501 : return RgX_Rg_mul(y2, gdiv(x,y1));
2280 : }
2281 : static GEN
2282 1186733 : div_rfrac(GEN x, GEN y)
2283 1186733 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2284 :
2285 : /* x != 0 */
2286 : static GEN
2287 1338593 : div_ser_scal(GEN x, GEN t)
2288 : {
2289 1338593 : if (ser_isexactzero(x))
2290 : {
2291 28 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gdiv(gel(x,2), t), varn(x), 1);
2292 28 : setvalser(z, valser(x)); return z;
2293 : }
2294 6202113 : pari_APPLY_ser(gdiv(gel(x,i), t));
2295 : }
2296 : GEN
2297 7658 : ser_normalize(GEN x)
2298 : {
2299 7658 : long i, lx = lg(x);
2300 : GEN c, z;
2301 7658 : if (lx == 2) return x;
2302 7658 : c = gel(x,2); if (gequal1(c)) return x;
2303 7581 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2304 108752 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2305 7581 : return z;
2306 : }
2307 :
2308 : /* y != 0 */
2309 : static GEN
2310 6694004 : div_T_scal(GEN x, GEN y, long tx) {
2311 6694004 : switch(tx)
2312 : {
2313 5344328 : case t_POL: return RgX_Rg_div(x, y);
2314 1338586 : case t_SER: return div_ser_scal(x, y);
2315 11090 : case t_RFRAC: return div_rfrac_scal(x,y);
2316 : }
2317 0 : pari_err_TYPE2("/",x,y);
2318 : return NULL; /* LCOV_EXCL_LINE */
2319 : }
2320 :
2321 : static GEN
2322 9263817 : div_scal_pol(GEN x, GEN y) {
2323 9263817 : long ly = lg(y);
2324 : pari_sp av;
2325 : GEN _1;
2326 9263817 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2327 9184737 : if (isrationalzero(x)) return zeropol(varn(y));
2328 7130641 : av = avma;
2329 7130641 : _1 = Rg_get_1(x); if (transtype(_1)) y = gmul(y, _1);
2330 7130641 : return gerepileupto(av, gred_rfrac_simple(x,y));
2331 : }
2332 : static GEN
2333 18550 : div_scal_ser(GEN x, GEN y)
2334 : {
2335 18550 : pari_sp av = avma;
2336 18550 : GEN _1 = Rg_get_1(x);
2337 18550 : if (transtype(_1)) y = gmul(y, _1);
2338 18550 : return gerepileupto(av, gmul(x, ser_inv(y)));
2339 : }
2340 : static GEN
2341 9271983 : div_scal_T(GEN x, GEN y, long ty) {
2342 9271983 : switch(ty)
2343 : {
2344 9216681 : case t_POL: return div_scal_pol(x, y);
2345 18543 : case t_SER: return div_scal_ser(x, y);
2346 36759 : case t_RFRAC: return div_scal_rfrac(x, y);
2347 : }
2348 0 : pari_err_TYPE2("/",x,y);
2349 : return NULL; /* LCOV_EXCL_LINE */
2350 : }
2351 :
2352 : /* assume tx = ty = t_SER, same variable vx */
2353 : static GEN
2354 771510 : div_ser(GEN x, GEN y, long vx)
2355 : {
2356 771510 : long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
2357 771510 : GEN y0 = y, z;
2358 771510 : pari_sp av = avma;
2359 :
2360 771510 : if (!signe(y)) pari_err_INV("div_ser", y);
2361 771503 : if (ser_isexactzero(x))
2362 : {
2363 59892 : if (lx == 2) return zeroser(vx, v);
2364 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2365 147 : setvalser(z, v); return z;
2366 : }
2367 711611 : if (lx < ly) ly = lx;
2368 711611 : y = ser2pol_i_normalize(y, ly, &e);
2369 711611 : if (e)
2370 : {
2371 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2372 7 : ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
2373 : }
2374 711611 : z = init_ser(ly, vx, v);
2375 711611 : if (ly == 3)
2376 : {
2377 12347 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2378 12347 : if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
2379 12347 : return gerepileupto(av, z);
2380 : }
2381 699264 : x = ser2pol_i(x, ly);
2382 699264 : y = RgXn_div_i(x, y, ly-2);
2383 699264 : return gc_GEN(av, fill_ser(z,y));
2384 : }
2385 : /* x,y compatible PADIC */
2386 : static GEN
2387 13212186 : divpp(GEN x, GEN y)
2388 : {
2389 13212186 : GEN M, ux = padic_u(x), uy = padic_u(y), p = padic_p(x);
2390 13212186 : long a = precp(x), b = precp(y), v = valp(x) - valp(y);
2391 :
2392 13212186 : if (!signe(uy)) pari_err_INV("divpp",y);
2393 13212183 : if (!signe(ux)) return zeropadic(p, v);
2394 13211854 : if (a > b) M = padic_pd(y); else { M = padic_pd(x); b = a; }
2395 13211854 : retmkpadic(Fp_div(ux, uy, M), icopy(padic_p(x)), icopy(M), v, b);
2396 : }
2397 : static GEN
2398 37093 : div_polmod_same(GEN T, GEN x, GEN y)
2399 : {
2400 37093 : long v = varn(T);
2401 37093 : GEN a, z = cgetg(3, t_POLMOD);
2402 37093 : gel(z,1) = RgX_copy(T);
2403 37093 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2404 11417 : a = gdiv(x, y);
2405 25676 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2406 119 : {
2407 119 : pari_sp av = avma;
2408 119 : a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
2409 : }
2410 25557 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2411 5278 : {
2412 5278 : pari_sp av = avma;
2413 5278 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2414 5278 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2415 5278 : a = gerepileupto(av, a);
2416 : }
2417 : else
2418 : {
2419 20279 : pari_sp av = avma;
2420 20279 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2421 20279 : a = gerepileupto(av, a);
2422 : }
2423 37093 : gel(z,2) = a; return z;
2424 : }
2425 : GEN
2426 432935317 : gdiv(GEN x, GEN y)
2427 : {
2428 432935317 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2429 : pari_sp av, tetpil;
2430 : GEN z, p1, p2;
2431 :
2432 432935317 : if (tx == ty) switch(tx)
2433 : {
2434 88718917 : case t_INT:
2435 88718917 : return Qdivii(x,y);
2436 :
2437 93660907 : case t_REAL: return divrr(x,y);
2438 25755 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2439 25755 : z = cgetg(3,t_INTMOD);
2440 25755 : if (X==Y || equalii(X,Y))
2441 25741 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2442 14 : gel(z,1) = gcdii(X,Y);
2443 14 : warn_coercion(X,Y,gel(z,1));
2444 14 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2445 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
2446 : }
2447 4041746 : case t_FRAC: {
2448 4041746 : GEN x1 = gel(x,1), x2 = gel(x,2);
2449 4041746 : GEN y1 = gel(y,1), y2 = gel(y,2);
2450 4041746 : z = cgetg(3, t_FRAC);
2451 4041746 : p1 = gcdii(x1, y1);
2452 4041744 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2453 4041744 : p1 = gcdii(x2, y2);
2454 4041745 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2455 4041744 : tetpil = avma;
2456 4041744 : gel(z,2) = mulii(x2,y1);
2457 4041743 : gel(z,1) = mulii(x1,y2);
2458 4041745 : normalize_frac(z);
2459 4041745 : fix_frac_if_int_GC(z,tetpil);
2460 4041746 : return z;
2461 : }
2462 9201112 : case t_COMPLEX:
2463 9201112 : if (isintzero(gel(y,1)))
2464 : {
2465 195412 : y = gel(y,2);
2466 195412 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2467 116011 : z = cgetg(3,t_COMPLEX);
2468 116011 : gel(z,1) = gdiv(gel(x,2), y);
2469 116011 : av = avma;
2470 116011 : gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
2471 116011 : return z;
2472 : }
2473 9005698 : av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
2474 9005721 : return gerepile(av, tetpil, gdiv(p2,p1));
2475 :
2476 587462 : case t_PADIC:
2477 587462 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("/",x,y);
2478 587462 : return divpp(x, y);
2479 :
2480 238 : case t_QUAD:
2481 238 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2482 224 : av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
2483 224 : return gerepile(av, tetpil, gdiv(p2,p1));
2484 :
2485 245921 : case t_FFELT: return FF_div(x,y);
2486 :
2487 41468 : case t_POLMOD:
2488 41468 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2489 37093 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2490 : else {
2491 4375 : av = avma;
2492 4375 : z = gerepileupto(av, gmul(x, ginv(y)));
2493 : }
2494 41468 : return z;
2495 :
2496 18549494 : case t_POL:
2497 18549494 : vx = varn(x);
2498 18549494 : vy = varn(y);
2499 18549494 : if (vx != vy) {
2500 154994 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2501 47136 : else return div_scal_pol(x, y);
2502 : }
2503 18394500 : if (!signe(y)) pari_err_INV("gdiv",y);
2504 18394500 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2505 18266545 : av = avma;
2506 18266545 : return gerepileupto(av, gred_rfrac2(x,y));
2507 :
2508 771531 : case t_SER:
2509 771531 : vx = varn(x);
2510 771531 : vy = varn(y);
2511 771531 : if (vx != vy) {
2512 21 : if (varncmp(vx, vy) < 0)
2513 : {
2514 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2515 7 : return div_ser_scal(x, y);
2516 : }
2517 7 : return div_scal_ser(x, y);
2518 : }
2519 771510 : return div_ser(x, y, vx);
2520 1189421 : case t_RFRAC:
2521 1189421 : vx = varn(gel(x,2));
2522 1189421 : vy = varn(gel(y,2));
2523 1189421 : if (vx != vy) {
2524 2688 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2525 756 : else return div_scal_rfrac(x, y);
2526 : }
2527 1186733 : return div_rfrac(x,y);
2528 :
2529 21 : case t_VEC: /* handle extended t_QFB */
2530 21 : case t_QFB: av = avma; return gerepileupto(av, qfbcomp(x, ginv(y)));
2531 :
2532 28 : case t_MAT:
2533 28 : p1 = RgM_div(x,y);
2534 28 : if (!p1) pari_err_INV("gdiv",y);
2535 21 : return p1;
2536 :
2537 0 : default: pari_err_TYPE2("/",x,y);
2538 : }
2539 :
2540 215906110 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2541 : {
2542 3764808 : long s = signe(x);
2543 3764808 : if (!s) {
2544 514660 : if (gequal0(y)) pari_err_INV("gdiv",y);
2545 514660 : switch (ty)
2546 : {
2547 511496 : default: return gen_0;
2548 28 : case t_INTMOD:
2549 28 : z = cgetg(3,t_INTMOD);
2550 28 : gel(z,1) = icopy(gel(y,1));
2551 28 : gel(z,2) = gen_0; return z;
2552 3136 : case t_FFELT: return FF_zero(y);
2553 : }
2554 : }
2555 3250148 : if (is_pm1(x)) {
2556 1302513 : if (s > 0) return ginv(y);
2557 221342 : av = avma; return gerepileupto(av, ginv(gneg(y)));
2558 : }
2559 1947635 : switch(ty)
2560 : {
2561 639572 : case t_REAL: return divir(x,y);
2562 21 : case t_INTMOD:
2563 21 : z = cgetg(3, t_INTMOD);
2564 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2565 792270 : case t_FRAC:
2566 792270 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2567 792270 : if (equali1(p1))
2568 : {
2569 271449 : set_avma((pari_sp)z);
2570 271449 : gel(z,2) = icopy(gel(y,1));
2571 271449 : gel(z,1) = mulii(gel(y,2), x);
2572 271449 : normalize_frac(z);
2573 271449 : fix_frac_if_int(z);
2574 : }
2575 : else
2576 : {
2577 520821 : x = diviiexact(x,p1); tetpil = avma;
2578 520821 : gel(z,2) = diviiexact(gel(y,1), p1);
2579 520821 : gel(z,1) = mulii(gel(y,2), x);
2580 520821 : normalize_frac(z);
2581 520821 : fix_frac_if_int_GC(z,tetpil);
2582 : }
2583 792270 : return z;
2584 :
2585 462 : case t_FFELT: return Z_FF_div(x,y);
2586 513579 : case t_COMPLEX: return divRc(x,y);
2587 1505 : case t_PADIC: return divTp(x, y);
2588 231 : case t_QUAD:
2589 231 : av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
2590 231 : return gerepile(av, tetpil, gdiv(p2,p1));
2591 : }
2592 : }
2593 212141298 : if (gequal0(y))
2594 : {
2595 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2596 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2597 : }
2598 :
2599 212151014 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2600 : {
2601 172826261 : case t_REAL:
2602 172826261 : switch(ty)
2603 : {
2604 170460135 : case t_INT: return divri(x,y);
2605 522601 : case t_FRAC:
2606 522601 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2607 522601 : return gerepileuptoleaf(av, z);
2608 1843589 : case t_COMPLEX: return divRc(x, y);
2609 42 : case t_QUAD: return divfq(x, y, realprec(x));
2610 18 : default: pari_err_TYPE2("/",x,y);
2611 : }
2612 :
2613 280 : case t_INTMOD:
2614 : switch(ty)
2615 : {
2616 196 : case t_INT:
2617 196 : z = cgetg(3, t_INTMOD);
2618 196 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2619 28 : case t_FRAC: { GEN X = gel(x,1);
2620 28 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2621 28 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2622 : }
2623 7 : case t_FFELT:
2624 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2625 0 : pari_err_OP("/",x,y);
2626 7 : return Z_FF_div(gel(x,2),y);
2627 :
2628 0 : case t_COMPLEX:
2629 0 : av = avma;
2630 0 : return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
2631 :
2632 14 : case t_QUAD:
2633 14 : av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
2634 14 : return gerepile(av,tetpil, gdiv(p2,p1));
2635 :
2636 7 : case t_PADIC: { GEN X = gel(x,1);
2637 7 : z = cgetg(3, t_INTMOD);
2638 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2639 : }
2640 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2641 : }
2642 :
2643 : case t_FRAC:
2644 : switch(ty)
2645 : {
2646 2099590 : case t_INT: z = cgetg(3, t_FRAC);
2647 2099590 : p1 = gcdii(y,gel(x,1));
2648 2099590 : if (equali1(p1))
2649 : {
2650 830124 : set_avma((pari_sp)z); tetpil = 0;
2651 830124 : gel(z,1) = icopy(gel(x,1));
2652 : }
2653 : else
2654 : {
2655 1269466 : y = diviiexact(y,p1); tetpil = avma;
2656 1269466 : gel(z,1) = diviiexact(gel(x,1), p1);
2657 : }
2658 2099590 : gel(z,2) = mulii(gel(x,2),y);
2659 2099590 : normalize_frac(z);
2660 2099590 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2661 2099590 : return z;
2662 :
2663 59530 : case t_REAL:
2664 59530 : av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
2665 59530 : return gerepile(av, tetpil, divir(gel(x,1), p1));
2666 :
2667 7 : case t_INTMOD: { GEN Y = gel(y,1);
2668 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2669 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2670 : }
2671 :
2672 28 : case t_FFELT: av=avma;
2673 28 : return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2674 :
2675 868 : case t_COMPLEX: return divRc(x, y);
2676 :
2677 2141 : case t_PADIC:
2678 2141 : if (!signe(gel(x,1))) return gen_0;
2679 2141 : return divTp(x, y);
2680 :
2681 42 : case t_QUAD:
2682 42 : av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
2683 42 : return gerepile(av,tetpil,gdiv(p2,p1));
2684 : }
2685 :
2686 : case t_FFELT:
2687 161 : switch (ty)
2688 : {
2689 77 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2690 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2691 7 : case t_INTMOD:
2692 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2693 0 : pari_err_OP("/",x,y);
2694 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2695 49 : default:
2696 49 : pari_err_TYPE2("/",x,y);
2697 : }
2698 0 : break;
2699 :
2700 13513957 : case t_COMPLEX:
2701 : switch(ty)
2702 : {
2703 13513958 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2704 0 : case t_INTMOD: return mulRc_direct(ginv(y), x);
2705 0 : case t_PADIC:
2706 0 : return Zp_nosquare_m1(padic_p(y))? divcR(x,y): divTp(x, y);
2707 0 : case t_QUAD:
2708 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2709 0 : return divfq(x, y, lx);
2710 : }
2711 :
2712 : case t_PADIC:
2713 : switch(ty)
2714 : {
2715 1205001 : case t_INT: case t_FRAC: { GEN p = padic_p(x);
2716 1204910 : return signe(padic_u(x))? divpT(x, y)
2717 2409911 : : zeropadic(p, valp(x) - Q_pval(y,p));
2718 : }
2719 7 : case t_INTMOD: { GEN Y = gel(y,1);
2720 7 : z = cgetg(3, t_INTMOD);
2721 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2722 : }
2723 14 : case t_COMPLEX: case t_QUAD:
2724 14 : av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
2725 14 : return gerepile(av,tetpil,gdiv(p1,p2));
2726 :
2727 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2728 : }
2729 :
2730 : case t_QUAD:
2731 : switch (ty)
2732 : {
2733 1190 : case t_INT: case t_INTMOD: case t_FRAC:
2734 1190 : z = cgetg(4,t_QUAD);
2735 1190 : gel(z,1) = ZX_copy(gel(x,1));
2736 1190 : gel(z,2) = gdiv(gel(x,2), y);
2737 1190 : gel(z,3) = gdiv(gel(x,3), y); return z;
2738 63 : case t_REAL: return divqf(x, y, realprec(y));
2739 14 : case t_PADIC: return divTp(x, y);
2740 0 : case t_COMPLEX:
2741 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2742 0 : return divqf(x, y, ly);
2743 : }
2744 : }
2745 22441617 : switch(ty) {
2746 582279 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2747 582279 : return gmul(x, ginv(y)); /* missing gerepile, for speed */
2748 42 : case t_MAT:
2749 42 : av = avma; p1 = RgM_inv(y);
2750 35 : if (!p1) pari_err_INV("gdiv",y);
2751 28 : return gerepileupto(av, gmul(x, p1));
2752 0 : case t_VEC: case t_COL:
2753 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2754 0 : pari_err_TYPE2("/",x,y);
2755 : }
2756 21859296 : switch(tx) {
2757 3791918 : case t_VEC: case t_COL: case t_MAT:
2758 12819027 : pari_APPLY_same(gdiv(gel(x,i),y));
2759 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2760 0 : pari_err_TYPE2("/",x,y);
2761 : }
2762 :
2763 18067378 : vy = gvar(y);
2764 18068392 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2765 209777 : vx = varn(X);
2766 209777 : if (vx != vy) {
2767 209203 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2768 208790 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2769 : }
2770 : /* y is POL, SER or RFRAC */
2771 574 : av = avma;
2772 574 : switch(ty)
2773 : {
2774 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2775 574 : default: y = ginvmod(gmod(y,X), X);
2776 : }
2777 567 : return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
2778 : }
2779 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2780 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2781 : * variable NO_VARIABLE, then the operation is incorrect */
2782 17858615 : vx = gvar(x);
2783 17858610 : if (vx != vy) { /* includes cases where one is scalar */
2784 15965582 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2785 9271573 : else return div_scal_T(x, y, ty);
2786 : }
2787 1893028 : switch(tx)
2788 : {
2789 1046171 : case t_POL:
2790 : switch(ty)
2791 : {
2792 28 : case t_SER:
2793 : {
2794 28 : GEN y0 = y;
2795 : long v;
2796 28 : av = avma; v = RgX_valrem(x, &x);
2797 28 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2798 14 : v -= valser(y); ly = lg(y); /* > 2 */
2799 14 : y = ser2pol_i_normalize(y, ly, &i);
2800 14 : if (i)
2801 : {
2802 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2803 7 : ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
2804 : }
2805 7 : z = init_ser(ly, vx, v);
2806 7 : return gc_GEN(av, fill_ser(z, RgXn_div(x, y, ly-2)));
2807 : }
2808 :
2809 1046143 : case t_RFRAC:
2810 : {
2811 1046143 : GEN y1 = gel(y,1), y2 = gel(y,2);
2812 1046143 : if (typ(y1) == t_POL && varn(y1) == vx)
2813 1143 : return mul_rfrac_scal(y2, y1, x);
2814 1045000 : av = avma;
2815 1045000 : return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2816 : }
2817 : }
2818 0 : break;
2819 :
2820 584324 : case t_SER:
2821 : switch(ty)
2822 : {
2823 584310 : case t_POL:
2824 : {
2825 584310 : long v = valser(x);
2826 584310 : lx = lg(x);
2827 584310 : if (lx == 2) return zeroser(vx, v - RgX_val(y));
2828 584205 : av = avma;
2829 584205 : x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
2830 584205 : z = init_ser(lx, vx, v);
2831 584205 : if (!signe(x)) setsigne(z,0);
2832 584205 : return gc_GEN(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
2833 : }
2834 14 : case t_RFRAC:
2835 14 : av = avma;
2836 14 : return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2837 : }
2838 0 : break;
2839 :
2840 262512 : case t_RFRAC:
2841 : switch(ty)
2842 : {
2843 262512 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2844 0 : case t_SER:
2845 0 : av = avma;
2846 0 : return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2847 : }
2848 0 : break;
2849 : }
2850 21 : pari_err_TYPE2("/",x,y);
2851 : return NULL; /* LCOV_EXCL_LINE */
2852 : }
2853 :
2854 : /********************************************************************/
2855 : /** **/
2856 : /** SIMPLE MULTIPLICATION **/
2857 : /** **/
2858 : /********************************************************************/
2859 : GEN
2860 229686855 : gmulsg(long s, GEN x)
2861 : {
2862 : pari_sp av;
2863 : long i;
2864 : GEN z;
2865 :
2866 229686855 : switch(typ(x))
2867 : {
2868 136885283 : case t_INT: return mulsi(s,x);
2869 69473609 : case t_REAL: return s? mulsr(s,x): gen_0; /* gmul semantic */
2870 174387 : case t_INTMOD: { GEN p = gel(x,1);
2871 174387 : z = cgetg(3,t_INTMOD);
2872 174387 : gel(z,2) = gc_INT((pari_sp)z, modii(mulsi(s,gel(x,2)), p));
2873 174387 : gel(z,1) = icopy(p); return z;
2874 : }
2875 548098 : case t_FFELT: return FF_Z_mul(x,stoi(s));
2876 6210706 : case t_FRAC:
2877 6210706 : if (!s) return gen_0;
2878 6196006 : z = cgetg(3,t_FRAC);
2879 6196008 : i = labs(s); i = ugcdiu(gel(x,2), i);
2880 6196010 : if (i == 1)
2881 : {
2882 2269441 : gel(z,2) = icopy(gel(x,2));
2883 2269441 : gel(z,1) = mulis(gel(x,1), s);
2884 : }
2885 : else
2886 : {
2887 3926569 : gel(z,2) = diviuexact(gel(x,2), (ulong)i);
2888 3926569 : gel(z,1) = mulis(gel(x,1), s/i);
2889 3926569 : fix_frac_if_int(z);
2890 : }
2891 6196006 : return z;
2892 :
2893 14215214 : case t_COMPLEX:
2894 14215214 : if (!s) return gen_0;
2895 13089410 : z = cgetg(3, t_COMPLEX);
2896 13089353 : gel(z,1) = gmulsg(s,gel(x,1));
2897 13087907 : gel(z,2) = gmulsg(s,gel(x,2)); return z;
2898 :
2899 1449 : case t_PADIC:
2900 1449 : if (!s) return gen_0;
2901 1449 : av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),x), x));
2902 :
2903 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2904 7 : gel(z,1) = ZX_copy(gel(x,1));
2905 7 : gel(z,2) = gmulsg(s,gel(x,2));
2906 7 : gel(z,3) = gmulsg(s,gel(x,3)); return z;
2907 :
2908 204442 : case t_POLMOD:
2909 204442 : retmkpolmod(gmulsg(s,gel(x,2)), RgX_copy(gel(x,1)));
2910 :
2911 748905 : case t_POL:
2912 748905 : if (!signe(x)) return RgX_copy(x);
2913 732105 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2914 3096258 : pari_APPLY_pol(gmulsg(s,gel(x,i)));
2915 :
2916 182 : case t_SER:
2917 182 : if (ser_isexactzero(x)) return gcopy(x);
2918 182 : if (!s) return Rg_get_0(x);
2919 3864 : pari_APPLY_ser(gmulsg(s,gel(x,i)));
2920 :
2921 0 : case t_RFRAC:
2922 0 : if (!s) return zeropol(varn(gel(x,2)));
2923 0 : if (s == 1) return gcopy(x);
2924 0 : if (s == -1) return gneg(x);
2925 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), stoi(s));
2926 :
2927 1233828 : case t_VEC: case t_COL: case t_MAT:
2928 3897151 : pari_APPLY_same(gmulsg(s,gel(x,i)));
2929 : }
2930 0 : pari_err_TYPE("gmulsg",x);
2931 : return NULL; /* LCOV_EXCL_LINE */
2932 : }
2933 :
2934 : GEN
2935 124340561 : gmulug(ulong s, GEN x)
2936 : {
2937 : pari_sp av;
2938 : long i;
2939 : GEN z;
2940 :
2941 124340561 : switch(typ(x))
2942 : {
2943 122223429 : case t_INT: return mului(s,x);
2944 1059371 : case t_REAL: return s? mulur(s,x): gen_0; /* gmul semantic */
2945 364 : case t_INTMOD: { GEN p = gel(x,1);
2946 364 : z = cgetg(3,t_INTMOD);
2947 364 : gel(z,2) = gc_INT((pari_sp)z, modii(mului(s,gel(x,2)), p));
2948 364 : gel(z,1) = icopy(p); return z;
2949 : }
2950 413 : case t_FFELT: return FF_Z_mul(x,utoi(s));
2951 981585 : case t_FRAC:
2952 981585 : if (!s) return gen_0;
2953 981571 : z = cgetg(3,t_FRAC);
2954 981571 : i = ugcdiu(gel(x,2), s);
2955 981571 : if (i == 1)
2956 : {
2957 822559 : gel(z,2) = icopy(gel(x,2));
2958 822559 : gel(z,1) = muliu(gel(x,1), s);
2959 : }
2960 : else
2961 : {
2962 159012 : gel(z,2) = diviuexact(gel(x,2), i);
2963 159012 : gel(z,1) = muliu(gel(x,1), s/i);
2964 159012 : fix_frac_if_int(z);
2965 : }
2966 981571 : return z;
2967 :
2968 30765 : case t_COMPLEX:
2969 30765 : if (!s) return gen_0;
2970 30765 : z = cgetg(3, t_COMPLEX);
2971 30765 : gel(z,1) = gmulug(s,gel(x,1));
2972 30765 : gel(z,2) = gmulug(s,gel(x,2)); return z;
2973 :
2974 7342 : case t_PADIC:
2975 7342 : if (!s) return gen_0;
2976 7342 : av = avma; return gerepileupto(av, mulpp(cvtop2(utoi(s),x), x));
2977 :
2978 0 : case t_QUAD: z = cgetg(4, t_QUAD);
2979 0 : gel(z,1) = ZX_copy(gel(x,1));
2980 0 : gel(z,2) = gmulug(s,gel(x,2));
2981 0 : gel(z,3) = gmulug(s,gel(x,3)); return z;
2982 :
2983 6783 : case t_POLMOD:
2984 6783 : retmkpolmod(gmulug(s,gel(x,2)), RgX_copy(gel(x,1)));
2985 :
2986 16226 : case t_POL:
2987 16226 : if (!signe(x)) return RgX_copy(x);
2988 15449 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2989 52843 : pari_APPLY_pol(gmulug(s,gel(x,i)));
2990 :
2991 0 : case t_SER:
2992 0 : if (ser_isexactzero(x)) return gcopy(x);
2993 0 : if (!s) return Rg_get_0(x);
2994 0 : pari_APPLY_ser(gmulug(s,gel(x,i)));
2995 :
2996 0 : case t_RFRAC:
2997 0 : if (!s) return zeropol(varn(gel(x,2)));
2998 0 : if (s == 1) return gcopy(x);
2999 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), utoi(s));
3000 :
3001 14286 : case t_VEC: case t_COL: case t_MAT:
3002 74458 : pari_APPLY_same(gmulug(s,gel(x,i)));
3003 : }
3004 0 : pari_err_TYPE("gmulsg",x);
3005 : return NULL; /* LCOV_EXCL_LINE */
3006 : }
3007 :
3008 : /********************************************************************/
3009 : /** **/
3010 : /** SIMPLE DIVISION **/
3011 : /** **/
3012 : /********************************************************************/
3013 :
3014 : GEN
3015 12965965 : gdivgs(GEN x, long s)
3016 : {
3017 12965965 : long tx = typ(x), i;
3018 : pari_sp av;
3019 : GEN z;
3020 :
3021 12965965 : if (!s)
3022 : {
3023 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3024 0 : pari_err_INV("gdivgs",gen_0);
3025 : }
3026 12965967 : switch(tx)
3027 : {
3028 1583916 : case t_INT: return Qdivis(x, s);
3029 8503614 : case t_REAL: return divrs(x,s);
3030 :
3031 357 : case t_INTMOD:
3032 357 : z = cgetg(3, t_INTMOD);
3033 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3034 :
3035 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3036 :
3037 552267 : case t_FRAC: z = cgetg(3, t_FRAC);
3038 552267 : i = ugcdiu(gel(x,1), labs(s));
3039 552267 : if (i == 1)
3040 : {
3041 409425 : gel(z,2) = mulsi(s, gel(x,2));
3042 409425 : gel(z,1) = icopy(gel(x,1));
3043 : }
3044 : else
3045 : {
3046 142842 : gel(z,2) = mulsi(s/i, gel(x,2));
3047 142842 : gel(z,1) = divis(gel(x,1), i);
3048 : }
3049 552267 : normalize_frac(z);
3050 552267 : fix_frac_if_int(z); return z;
3051 :
3052 1841738 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3053 1841738 : gel(z,1) = gdivgs(gel(x,1),s);
3054 1841737 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3055 :
3056 133 : case t_PADIC: /* divpT */
3057 : {
3058 133 : GEN p = padic_p(x);
3059 133 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3060 133 : av = avma;
3061 133 : return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
3062 : }
3063 :
3064 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3065 28 : gel(z,1) = ZX_copy(gel(x,1));
3066 28 : gel(z,2) = gdivgs(gel(x,2),s);
3067 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3068 :
3069 37849 : case t_POLMOD:
3070 37849 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3071 :
3072 91 : case t_RFRAC:
3073 91 : if (s == 1) return gcopy(x);
3074 84 : else if (s == -1) return gneg(x);
3075 84 : return div_rfrac_scal(x, stoi(s));
3076 :
3077 157596 : case t_POL: pari_APPLY_pol_normalized(gdivgs(gel(x,i),s));
3078 0 : case t_SER: pari_APPLY_ser_normalized(gdivgs(gel(x,i),s));
3079 407593 : case t_VEC:
3080 : case t_COL:
3081 894977 : case t_MAT: pari_APPLY_same(gdivgs(gel(x,i),s));
3082 : }
3083 0 : pari_err_TYPE2("/",x, stoi(s));
3084 : return NULL; /* LCOV_EXCL_LINE */
3085 : }
3086 :
3087 : GEN
3088 54716433 : gdivgu(GEN x, ulong s)
3089 : {
3090 54716433 : long tx = typ(x), i;
3091 : pari_sp av;
3092 : GEN z;
3093 :
3094 54716433 : if (!s)
3095 : {
3096 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3097 0 : pari_err_INV("gdivgu",gen_0);
3098 : }
3099 54716443 : switch(tx)
3100 : {
3101 17916052 : case t_INT: return Qdiviu(x, s);
3102 12957167 : case t_REAL: return divru(x,s);
3103 :
3104 210315 : case t_INTMOD:
3105 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3106 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3107 :
3108 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3109 :
3110 642010 : case t_FRAC: z = cgetg(3, t_FRAC);
3111 642010 : i = ugcdiu(gel(x,1), s);
3112 642010 : if (i == 1)
3113 : {
3114 462134 : gel(z,2) = mului(s, gel(x,2));
3115 462134 : gel(z,1) = icopy(gel(x,1));
3116 : }
3117 : else
3118 : {
3119 179876 : gel(z,2) = mului(s/i, gel(x,2));
3120 179876 : gel(z,1) = divis(gel(x,1), i);
3121 : }
3122 642010 : normalize_frac(z);
3123 642010 : fix_frac_if_int(z); return z;
3124 :
3125 11384779 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3126 11384779 : gel(z,1) = gdivgu(gel(x,1),s);
3127 11384778 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3128 :
3129 11551079 : case t_PADIC: /* divpT */
3130 : {
3131 11551079 : GEN p = padic_p(x);
3132 11551079 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3133 11415487 : av = avma;
3134 11415487 : return gerepileupto(av, divpp(x, cvtop2(utoi(s),x)));
3135 : }
3136 :
3137 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3138 0 : gel(z,1) = ZX_copy(gel(x,1));
3139 0 : gel(z,2) = gdivgu(gel(x,2),s);
3140 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3141 :
3142 1456 : case t_POLMOD:
3143 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3144 :
3145 56 : case t_RFRAC:
3146 56 : if (s == 1) return gcopy(x);
3147 56 : return div_rfrac_scal(x, utoi(s));
3148 :
3149 205933 : case t_POL: pari_APPLY_pol_normalized(gdivgu(gel(x,i),s));
3150 16730 : case t_SER: pari_APPLY_ser_normalized(gdivgu(gel(x,i),s));
3151 336 : case t_VEC:
3152 : case t_COL:
3153 1162 : case t_MAT: pari_APPLY_same(gdivgu(gel(x,i),s));
3154 : }
3155 0 : pari_err_TYPE2("/",x, utoi(s));
3156 : return NULL; /* LCOV_EXCL_LINE */
3157 : }
3158 :
3159 : /* x / (i*(i+1)) */
3160 : GEN
3161 224299268 : divrunextu(GEN x, ulong i)
3162 : {
3163 224299268 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3164 0 : return divri(x, muluu(i , i+1));
3165 : else
3166 224299268 : return divru(x, i*(i+1));
3167 : }
3168 : /* x / (i*(i+1)) */
3169 : GEN
3170 1036796 : gdivgunextu(GEN x, ulong i)
3171 : {
3172 1036796 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3173 0 : return gdivgu(x, i*(i+1));
3174 : else
3175 1036796 : return gdiv(x, muluu(i, i+1));
3176 : }
3177 :
3178 : /* True shift (exact multiplication by 2^n) */
3179 : GEN
3180 133214499 : gmul2n(GEN x, long n)
3181 : {
3182 : GEN z, a, b;
3183 : long k, l;
3184 :
3185 133214499 : switch(typ(x))
3186 : {
3187 36353099 : case t_INT:
3188 36353099 : if (n>=0) return shifti(x,n);
3189 5049726 : if (!signe(x)) return gen_0;
3190 3353381 : l = vali(x); n = -n;
3191 3353408 : if (n<=l) return shifti(x,-n);
3192 398507 : z = cgetg(3,t_FRAC);
3193 398507 : gel(z,1) = shifti(x,-l);
3194 398507 : gel(z,2) = int2n(n-l); return z;
3195 :
3196 62831833 : case t_REAL:
3197 62831833 : return shiftr(x,n);
3198 :
3199 181013 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3200 181013 : z = cgetg(3,t_INTMOD);
3201 181013 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3202 82852 : gel(z,2) = gc_INT((pari_sp)z, modii(shifti(a,n), b));
3203 82850 : gel(z,1) = icopy(b); return z;
3204 :
3205 217424 : case t_FFELT: return FF_mul2n(x,n);
3206 :
3207 1271345 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3208 1271345 : l = vali(a);
3209 1271345 : k = vali(b);
3210 1271345 : if (n+l >= k)
3211 : {
3212 396730 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3213 303941 : l = n-k; k = -k;
3214 : }
3215 : else
3216 : {
3217 874615 : k = -(l+n); l = -l;
3218 : }
3219 1178556 : z = cgetg(3,t_FRAC);
3220 1178556 : gel(z,1) = shifti(a,l);
3221 1178556 : gel(z,2) = shifti(b,k); return z;
3222 :
3223 10728177 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3224 10728167 : gel(z,1) = gmul2n(gel(x,1),n);
3225 10728164 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3226 :
3227 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3228 105 : gel(z,1) = ZX_copy(gel(x,1));
3229 105 : gel(z,2) = gmul2n(gel(x,2),n);
3230 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3231 :
3232 193830 : case t_POLMOD:
3233 193830 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3234 :
3235 5063369 : case t_POL:
3236 21511722 : pari_APPLY_pol(gmul2n(gel(x,i),n));
3237 102847 : case t_SER:
3238 102847 : if (ser_isexactzero(x)) return gcopy(x);
3239 634256 : pari_APPLY_ser(gmul2n(gel(x,i),n));
3240 16267624 : case t_VEC: case t_COL: case t_MAT:
3241 59944623 : pari_APPLY_same(gmul2n(gel(x,i),n));
3242 :
3243 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3244 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3245 :
3246 5628 : case t_PADIC: /* int2n wrong if n < 0 */
3247 5628 : return gmul(gmul2n(gen_1,n),x);
3248 : }
3249 0 : pari_err_TYPE("gmul2n",x);
3250 : return NULL; /* LCOV_EXCL_LINE */
3251 : }
3252 :
3253 : /*******************************************************************/
3254 : /* */
3255 : /* INVERSE */
3256 : /* */
3257 : /*******************************************************************/
3258 : static GEN
3259 209968 : inv_polmod(GEN T, GEN x)
3260 : {
3261 209968 : GEN z = cgetg(3,t_POLMOD), a;
3262 209968 : gel(z,1) = RgX_copy(T);
3263 209967 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3264 85327 : a = ginv(x);
3265 : else
3266 : {
3267 124640 : if (lg(T) == 5) /* quadratic fields */
3268 13062 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3269 : else
3270 111578 : a = RgXQ_inv(x, T);
3271 : }
3272 209967 : gel(z,2) = a; return z;
3273 : }
3274 :
3275 : GEN
3276 36372218 : ginv(GEN x)
3277 : {
3278 : long s;
3279 : pari_sp av, tetpil;
3280 : GEN z, y, p1, p2;
3281 :
3282 36372218 : switch(typ(x))
3283 : {
3284 9404300 : case t_INT:
3285 9404300 : if (is_pm1(x)) return icopy(x);
3286 7883218 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3287 7883204 : z = cgetg(3,t_FRAC);
3288 7889696 : gel(z,1) = s<0? gen_m1: gen_1;
3289 7889696 : gel(z,2) = absi(x); return z;
3290 :
3291 10745114 : case t_REAL: return invr(x);
3292 :
3293 6062 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3294 6062 : gel(z,1) = icopy(gel(x,1));
3295 6062 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3296 :
3297 532072 : case t_FRAC: {
3298 532072 : GEN a = gel(x,1), b = gel(x,2);
3299 532072 : s = signe(a);
3300 532072 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3301 247533 : z = cgetg(3,t_FRAC);
3302 247534 : gel(z,1) = icopy(b);
3303 247534 : gel(z,2) = icopy(a);
3304 247534 : normalize_frac(z); return z;
3305 : }
3306 9614282 : case t_COMPLEX:
3307 9614282 : av=avma;
3308 9614282 : p1=cxnorm(x);
3309 9613482 : p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
3310 9614235 : tetpil=avma;
3311 9614235 : return gerepile(av,tetpil,divcR(p2,p1));
3312 :
3313 273 : case t_QUAD:
3314 273 : av=avma; p1=quadnorm(x); p2=conj_i(x); tetpil=avma;
3315 273 : return gerepile(av,tetpil,gdiv(p2,p1));
3316 :
3317 358358 : case t_PADIC:
3318 : {
3319 358358 : GEN p = padic_p(x), pd = padic_pd(x), u = padic_u(x);
3320 358358 : long d = precp(x);
3321 358358 : if (!signe(u)) pari_err_INV("ginv",x);
3322 358351 : retmkpadic(Zp_inv(u, p, d), icopy(p), icopy(pd), -valp(x), d);
3323 : }
3324 :
3325 209968 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3326 14031 : case t_FFELT: return FF_inv(x);
3327 5322891 : case t_POL: return gred_rfrac_simple(gen_1,x);
3328 26341 : case t_SER: return ser_inv(x);
3329 2933 : case t_RFRAC:
3330 : {
3331 2933 : GEN n = gel(x,1), d = gel(x,2);
3332 2933 : pari_sp av = avma, ltop;
3333 2933 : if (gequal0(n)) pari_err_INV("ginv",x);
3334 :
3335 2933 : n = simplify_shallow(n);
3336 2933 : if (typ(n) != t_POL || varn(n) != varn(d))
3337 : {
3338 2933 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3339 679 : ltop = avma;
3340 679 : z = RgX_Rg_div(d,n);
3341 : } else {
3342 0 : ltop = avma;
3343 0 : z = cgetg(3,t_RFRAC);
3344 0 : gel(z,1) = RgX_copy(d);
3345 0 : gel(z,2) = RgX_copy(n);
3346 : }
3347 679 : stackdummy(av, ltop);
3348 679 : return z;
3349 : }
3350 :
3351 21 : case t_VEC: if (!is_ext_qfr(x)) break;
3352 : case t_QFB:
3353 98338 : return qfbpow(x, gen_m1);
3354 38834 : case t_MAT:
3355 38834 : y = RgM_inv(x);
3356 38827 : if (!y) pari_err_INV("ginv",x);
3357 38757 : return y;
3358 28 : case t_VECSMALL:
3359 : {
3360 28 : long i, lx = lg(x)-1;
3361 28 : y = zero_zv(lx);
3362 112 : for (i=1; i<=lx; i++)
3363 : {
3364 84 : long xi = x[i];
3365 84 : if (xi<1 || xi>lx || y[xi])
3366 0 : pari_err_TYPE("ginv [not a permutation]", x);
3367 84 : y[xi] = i;
3368 : }
3369 28 : return y;
3370 : }
3371 : }
3372 6 : pari_err_TYPE("inverse",x);
3373 : return NULL; /* LCOV_EXCL_LINE */
3374 : }
|
