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 = gc_upto((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 = gc_upto((pari_sp)(z+3), gel(z,1));\
33 : else\
34 : gc_slice_unsafe((pari_sp)z, tetpil, z+1, 2); }
35 :
36 : static void
37 112 : warn_coercion(GEN x, GEN y, GEN z)
38 : {
39 112 : if (DEBUGLEVEL)
40 56 : pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
41 112 : }
42 :
43 : static long
44 35 : kro_quad(GEN x, GEN y)
45 35 : { 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 307902861 : addRc(GEN x, GEN y) {
61 307902861 : GEN z = cgetg(3,t_COMPLEX);
62 307884682 : gel(z,1) = gadd(x,gel(y,1));
63 307861967 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 354541989 : mulRc(GEN x, GEN y)
67 : {
68 354541989 : GEN a = gel(y,1), b = gel(y,2), z = cgetg(3,t_COMPLEX);
69 354430839 : gel(z,1) = (typ(x) != t_INTMOD && isintzero(a))? gen_0: gmul(x, a);
70 355115449 : gel(z,2) = gmul(x, b); return z;
71 : }
72 : static GEN
73 2286408 : divRc(GEN x, GEN y)
74 : {
75 2286408 : GEN a = gel(y,1), b = gel(y,2), z = cgetg(3,t_COMPLEX);
76 2286411 : pari_sp av = avma, av2;
77 2286411 : if (isintzero(a) && typ(x) != t_INTMOD)
78 : {
79 49294 : gel(z,1) = gen_0;
80 49294 : gel(z,2) = gc_upto(av, gdiv(x, gneg(b)));
81 : }
82 : else
83 : {
84 2237116 : GEN t = gdiv(x, cxnorm(y)), mb = gneg(b);
85 2237108 : av2 = avma;
86 2237108 : gel(z,1) = gmul(t, a);
87 2237100 : gel(z,2) = gmul(t, mb);
88 2237105 : gc_slice_unsafe(av, av2, z+1, 2);
89 : }
90 2286405 : return z;
91 : }
92 : static GEN
93 301 : divRq(GEN x, GEN y)
94 : {
95 301 : pari_sp av = avma;
96 301 : return gc_upto(av, gdiv(gmul(x,conj_i(y)), quadnorm(y)));
97 : }
98 : static GEN
99 23353235 : divcR(GEN x, GEN y) {
100 23353235 : GEN z = cgetg(3,t_COMPLEX);
101 23353178 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
102 23352952 : gel(z,2) = gdiv(gel(x,2), y); return z;
103 : }
104 : static GEN
105 1288 : addRq(GEN x, GEN y) {
106 1288 : GEN z = cgetg(4,t_QUAD);
107 1288 : gel(z,1) = ZX_copy(gel(y,1));
108 1288 : gel(z,2) = gadd(x, gel(y,2));
109 1288 : gel(z,3) = gcopy(gel(y,3)); return z;
110 : }
111 : static GEN
112 2324 : mulRq(GEN x, GEN y) {
113 2324 : GEN z = cgetg(4,t_QUAD);
114 2324 : gel(z,1) = ZX_copy(gel(y,1));
115 2324 : gel(z,2) = gmul(x,gel(y,2));
116 2324 : gel(z,3) = gmul(x,gel(y,3)); return z;
117 : }
118 : static GEN
119 77 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
120 77 : long i = gexpo(x) - gexpo(y);
121 77 : if (i > 0) prec += nbits2extraprec(i);
122 77 : return gc_upto(av, gadd(y, quadtofp(x, prec)));
123 : }
124 : static GEN
125 37741215 : mulrfrac(GEN x, GEN y)
126 : {
127 : pari_sp av;
128 37741215 : GEN z, a = gel(y,1), b = gel(y,2);
129 37741215 : if (is_pm1(a)) /* frequent special case */
130 : {
131 12722133 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
132 12722290 : return z;
133 : }
134 25018844 : av = avma;
135 25018844 : return gc_leaf(av, divri(mulri(x,a), b));
136 : }
137 : static GEN
138 28 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
139 28 : return gc_upto(av, gmul(y, quadtofp(x, prec)));
140 : }
141 : static GEN
142 63 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
143 63 : return gc_upto(av, gdiv(quadtofp(x,prec), y));
144 : }
145 : static GEN
146 42 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
147 42 : return gc_upto(av, gdiv(x, quadtofp(y,prec)));
148 : }
149 : /* y PADIC, x + y by converting x to padic */
150 : static GEN
151 7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
152 :
153 7 : if (!valp(y)) z = cvtop2(x,y);
154 : else {
155 7 : long l = signe(padic_u(y))? valp(y) + precp(y): valp(y);
156 7 : z = cvtop(x, padic_p(y), l);
157 : }
158 7 : return gc_upto(av, addsub_pp(z, y, addii));
159 : }
160 : /* y PADIC, x * y by converting x to padic */
161 : static GEN
162 4985917 : mulTp(GEN x, GEN y) { pari_sp av = avma;
163 4985917 : return gc_upto(av, mulpp(cvtop2(x,y), y));
164 : }
165 : /* y PADIC, non zero x / y by converting x to padic */
166 : static GEN
167 3674 : divTp(GEN x, GEN y) { pari_sp av = avma;
168 3674 : return gc_upto(av, divpp(cvtop2(x,y), y));
169 : }
170 : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
171 : * converted to O(p^e) and division by 0 */
172 : static GEN
173 1204910 : divpT(GEN x, GEN y) { pari_sp av = avma;
174 1204910 : return gc_upto(av, divpp(x, cvtop2(y,x)));
175 : }
176 :
177 : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
178 : * clean memory from z on */
179 : static GEN
180 3338891 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 3338891 : if (lgefint(X) == 3) {
182 3015761 : ulong u = Fl_add(itou(x),itou(y), X[2]);
183 3015761 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
184 : }
185 : else {
186 323130 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
187 323130 : gel(z,2) = gc_INT((pari_sp)z, u);
188 : }
189 3338891 : gel(z,1) = icopy(X); return z;
190 : }
191 : static GEN
192 1021769 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
193 1021769 : if (lgefint(X) == 3) {
194 647735 : ulong u = Fl_sub(itou(x),itou(y), X[2]);
195 647735 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
196 : }
197 : else {
198 374034 : GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
199 374034 : gel(z,2) = gc_INT((pari_sp)z, u);
200 : }
201 1021769 : gel(z,1) = icopy(X); return z;
202 : }
203 : /* cf add_intmod_same */
204 : static GEN
205 2788909 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
206 2788909 : if (lgefint(X) == 3) {
207 2647529 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
208 2647529 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 141380 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x,y), X) );
212 2788914 : gel(z,1) = icopy(X); return z;
213 : }
214 : /* cf add_intmod_same */
215 : static GEN
216 341509 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
217 : {
218 341509 : if (lgefint(X) == 3) {
219 318963 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
220 318956 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
221 : }
222 : else
223 22546 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
224 341503 : gel(z,1) = icopy(X); return z;
225 : }
226 :
227 : /*******************************************************************/
228 : /* */
229 : /* REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC) */
230 : /* */
231 : /* (static routines are not memory clean, but OK for gc_upto) */
232 : /*******************************************************************/
233 : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
234 : * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
235 : static GEN
236 9959022 : rfrac_denom_mul_scal(GEN d, GEN y)
237 : {
238 9959022 : GEN D = RgX_Rg_mul(d, y);
239 9959022 : if (lg(D) != lg(d))
240 : { /* try to generate a meaningful diagnostic */
241 0 : D = gdiv(leading_coeff(d), y); /* should fail */
242 0 : pari_err_INV("gred_rfrac", y); /* better than nothing */
243 : }
244 9959022 : return D;
245 : }
246 :
247 : static int
248 57790160 : Leading_is_neg(GEN x)
249 : {
250 122214139 : while (typ(x) == t_POL) x = leading_coeff(x);
251 57790160 : return is_real_t(typ(x))? gsigne(x) < 0: 0;
252 : }
253 :
254 : static int
255 154397901 : transtype(GEN x) { return x != gen_1 && typ(x) != t_PADIC; }
256 :
257 : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
258 : GEN
259 57794249 : gred_rfrac_simple(GEN n, GEN d)
260 : {
261 : GEN _1n, _1d, c, cn, cd, z;
262 57794249 : long dd = degpol(d);
263 :
264 57794249 : if (dd <= 0)
265 : {
266 4089 : if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
267 4089 : n = gdiv(n, gel(d,2));
268 4089 : if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
269 4089 : return n;
270 : }
271 57790160 : if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
272 57790160 : _1n = Rg_get_1(n);
273 57790160 : _1d = Rg_get_1(d);
274 57790160 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
275 57790160 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
276 57790153 : cd = content(d);
277 59659956 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
278 57790153 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
279 57790153 : if (!gequal1(cd)) {
280 6586881 : d = RgX_Rg_div(d,cd);
281 6586881 : if (!gequal1(cn))
282 : {
283 1313712 : if (gequal0(cn)) {
284 70 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
285 0 : n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
286 0 : c = gen_1;
287 : } else {
288 1313642 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
289 1313642 : c = gdiv(cn,cd);
290 : }
291 : }
292 : else
293 5273169 : c = ginv(cd);
294 : } else {
295 51203272 : if (!gequal1(cn))
296 : {
297 3341059 : if (gequal0(cn)) {
298 1008 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
299 21 : c = gen_1;
300 : } else {
301 3340051 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
302 3340051 : c = cn;
303 : }
304 : } else {
305 47862213 : GEN y = cgetg(3,t_RFRAC);
306 47862213 : gel(y,1) = gcopy(n);
307 47862213 : gel(y,2) = RgX_copy(d); return y;
308 : }
309 : }
310 :
311 9926883 : if (typ(c) == t_POL)
312 : {
313 912616 : z = c;
314 953342 : do { z = content(z); } while (typ(z) == t_POL);
315 912616 : cd = denom_i(z);
316 912616 : cn = gmul(c, cd);
317 : }
318 : else
319 : {
320 9014267 : cn = numer_i(c);
321 9014267 : cd = denom_i(c);
322 : }
323 9926883 : z = cgetg(3,t_RFRAC);
324 9926883 : gel(z,1) = gmul(n, cn);
325 9926883 : gel(z,2) = d = rfrac_denom_mul_scal(d, cd);
326 : /* can happen: Pol(O(17^-1)) / Pol([Mod(9,23), O(23^-3)]) */
327 9926883 : if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
328 9926883 : return z;
329 : }
330 :
331 : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
332 : static GEN
333 206780 : fix_rfrac(GEN x, long d)
334 : {
335 : GEN z, N, D;
336 206780 : if (!d || typ(x) == t_POL) return x;
337 165697 : z = cgetg(3, t_RFRAC);
338 165697 : N = gel(x,1);
339 165697 : D = gel(x,2);
340 165697 : if (d > 0) {
341 6650 : gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
342 160944 : : monomialcopy(N,d,varn(D));
343 160881 : gel(z, 2) = RgX_copy(D);
344 : } else {
345 4816 : gel(z, 1) = gcopy(N);
346 4816 : gel(z, 2) = RgX_shift(D, -d);
347 : }
348 165697 : return z;
349 : }
350 :
351 : /* assume d != 0 */
352 : static GEN
353 44823163 : gred_rfrac2(GEN n, GEN d)
354 : {
355 : GEN y, z, _1n, _1d;
356 : long v, vd, vn;
357 :
358 44823163 : n = simplify_shallow(n);
359 44823163 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
360 37965889 : d = simplify_shallow(d);
361 37965889 : if (typ(d) != t_POL) return gdiv(n,d);
362 36788993 : vd = varn(d);
363 36788993 : if (typ(n) != t_POL)
364 : {
365 20654572 : if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
366 20653165 : if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
367 0 : pari_err_BUG("gred_rfrac2 [incompatible variables]");
368 : }
369 16134421 : vn = varn(n);
370 16134421 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
371 15992470 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
372 15832697 : _1n = Rg_get_1(n);
373 15832697 : _1d = Rg_get_1(d);
374 15832697 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
375 15832697 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
376 :
377 : /* now n and d are t_POLs in the same variable */
378 15832697 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
379 15832697 : if (!degpol(d))
380 : {
381 12798208 : n = RgX_Rg_div(n,gel(d,2));
382 12798208 : return v? RgX_mulXn(n,v): n;
383 : }
384 :
385 : /* X does not divide gcd(n,d), deg(d) > 0 */
386 3034489 : if (!isinexact(n) && !isinexact(d))
387 : {
388 3034244 : y = RgX_divrem(n, d, &z);
389 3034244 : if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
390 206535 : z = RgX_gcd(d, z);
391 206535 : if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
392 : }
393 206780 : return fix_rfrac(gred_rfrac_simple(n,d), v);
394 : }
395 :
396 : /* x,y t_INT, return x/y in reduced form */
397 : GEN
398 141909138 : Qdivii(GEN x, GEN y)
399 : {
400 141909138 : pari_sp av = avma;
401 : GEN r, q;
402 :
403 141909138 : if (lgefint(y) == 3)
404 : {
405 123551990 : q = Qdiviu(x, y[2]);
406 123554269 : if (signe(y) > 0) return q;
407 11096510 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
408 11096538 : return q;
409 : }
410 18357148 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
411 18359448 : if (equali1(x))
412 : {
413 5186684 : if (!signe(y)) pari_err_INV("gdiv",y);
414 5186579 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
415 : }
416 13172785 : q = dvmdii(x,y,&r);
417 13172636 : if (r == gen_0) return q; /* gen_0 intended */
418 8735557 : r = gcdii(y, r);
419 8735603 : if (lgefint(r) == 3)
420 : {
421 7831153 : ulong t = r[2];
422 7831153 : set_avma(av);
423 7831131 : if (t == 1) q = mkfraccopy(x,y);
424 : else
425 : {
426 2831299 : q = cgetg(3,t_FRAC);
427 2831441 : gel(q,1) = diviuexact(x,t);
428 2831361 : gel(q,2) = diviuexact(y,t);
429 : }
430 : }
431 : else
432 : { /* rare: r and q left on stack for efficiency */
433 904450 : q = cgetg(3,t_FRAC);
434 904466 : gel(q,1) = diviiexact(x,r);
435 904467 : gel(q,2) = diviiexact(y,r);
436 : }
437 8735569 : normalize_frac(q); return q;
438 : }
439 :
440 : static GEN
441 64650224 : utoi_sign(ulong x, long s) { return s > 0? utoipos(x): utoineg(x); }
442 :
443 : /* x t_INT, return x/y in reduced form */
444 : GEN
445 148787897 : Qdiviu(GEN x, ulong y)
446 : {
447 : pari_sp av;
448 : ulong r, t;
449 : long s;
450 : GEN q;
451 :
452 148787897 : if (y == 1) return icopy(x);
453 126880108 : if (!y) pari_err_INV("gdiv",gen_0);
454 126885989 : s = signe(x); if (!s) return gen_0;
455 80996636 : if (lgefint(x) == 3)
456 : {
457 75175629 : ulong xx = uel(x,2), g;
458 75175629 : if (xx == 1) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
459 64648854 : g = ugcd(xx, y);
460 64650396 : if (g != 1)
461 : {
462 43292760 : xx /= g; if (y == g) return utoi_sign(xx, s);
463 16958845 : y /= g;
464 : }
465 38316481 : retmkfrac(utoi_sign(xx, s), utoipos(y));
466 : }
467 5821007 : av = avma; q = absdiviu_rem(x,y,&r);
468 5821483 : if (!r)
469 : {
470 2006052 : if (s < 0) togglesign(q);
471 2006052 : return q;
472 : }
473 3815431 : t = ugcd(y, r); set_avma(av);
474 3815531 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
475 1422708 : retmkfrac(diviuexact(x,t), utoipos(y / t));
476 : }
477 :
478 : /* x t_INT, return x/y in reduced form */
479 : GEN
480 1592288 : Qdivis(GEN x, long y)
481 : {
482 : GEN q;
483 1592288 : if (y >= 0) return Qdiviu(x, y);
484 180901 : q = Qdiviu(x, -y);
485 180901 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
486 180901 : return q;
487 : }
488 :
489 : /*******************************************************************/
490 : /* */
491 : /* CONJUGATION */
492 : /* */
493 : /*******************************************************************/
494 : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
495 : static GEN
496 18403 : quad_polmod_conj(GEN x, GEN y)
497 : {
498 : GEN z, u, v, a, b;
499 18403 : if (typ(x) != t_POL) return gcopy(x);
500 18403 : if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
501 18403 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
502 18403 : b = gel(y,3); v = gel(x,2);
503 18403 : z = cgetg(4, t_POL); z[1] = x[1];
504 18403 : gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
505 18403 : gel(z,3) = gneg(u); return z;
506 : }
507 : static GEN
508 18403 : quad_polmod_norm(GEN x, GEN y)
509 : {
510 : GEN z, u, v, a, b, c;
511 18403 : if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
512 18403 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
513 18403 : b = gel(y,3); v = gel(x,2);
514 18403 : c = gel(y,2);
515 18403 : z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
516 18403 : if (!gequal1(a)) z = gdiv(z, a);
517 18403 : return gadd(z, gsqr(v));
518 : }
519 :
520 : GEN
521 31492428 : conj_i(GEN x)
522 : {
523 31492428 : switch(typ(x))
524 : {
525 7459598 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
526 7459598 : return x;
527 :
528 23868210 : case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
529 :
530 945 : case t_QUAD:
531 : {
532 945 : GEN y = cgetg(4,t_QUAD);
533 945 : gel(y,1) = gel(x,1);
534 945 : gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
535 945 : : gadd(gel(x,2), gel(x,3));
536 945 : gel(y,3) = gneg(gel(x,3)); return y;
537 : }
538 350 : case t_POL: pari_APPLY_pol_normalized(conj_i(gel(x,i)));
539 32347 : case t_SER: pari_APPLY_ser_normalized(conj_i(gel(x,i)));
540 :
541 153418 : case t_RFRAC:
542 : case t_VEC:
543 : case t_COL:
544 553447 : case t_MAT: pari_APPLY_same(conj_i(gel(x,i)));
545 :
546 0 : case t_POLMOD:
547 : {
548 0 : GEN X = gel(x,1);
549 0 : long d = degpol(X);
550 0 : if (d < 2) return x;
551 0 : if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
552 : }
553 : }
554 0 : pari_err_TYPE("gconj",x);
555 : return NULL; /* LCOV_EXCL_LINE */
556 : }
557 : GEN
558 390135 : gconj(GEN x)
559 390135 : { pari_sp av = avma; return gc_GEN(av, conj_i(x)); }
560 :
561 : GEN
562 84 : conjvec(GEN x,long prec)
563 : {
564 : long lx, s, i;
565 : GEN z;
566 :
567 84 : switch(typ(x))
568 : {
569 0 : case t_INT: case t_INTMOD: case t_FRAC:
570 0 : return mkcolcopy(x);
571 :
572 0 : case t_COMPLEX: case t_QUAD:
573 0 : z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
574 :
575 28 : case t_FFELT:
576 28 : return FF_conjvec(x);
577 :
578 0 : case t_VEC: case t_COL:
579 0 : lx = lg(x); z = cgetg(lx,t_MAT);
580 0 : if (lx == 1) return z;
581 0 : gel(z,1) = conjvec(gel(x,1),prec);
582 0 : s = lgcols(z);
583 0 : for (i=2; i<lx; i++)
584 : {
585 0 : gel(z,i) = conjvec(gel(x,i),prec);
586 0 : if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
587 : }
588 0 : break;
589 :
590 56 : case t_POLMOD: {
591 56 : GEN T = gel(x,1), r;
592 : pari_sp av;
593 :
594 56 : lx = lg(T);
595 56 : if (lx <= 3) return cgetg(1,t_COL);
596 56 : x = gel(x,2);
597 238 : for (i=2; i<lx; i++)
598 : {
599 189 : GEN c = gel(T,i);
600 189 : switch(typ(c)) {
601 7 : case t_INTMOD: {
602 7 : GEN p = gel(c,1);
603 : pari_sp av;
604 7 : if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
605 7 : av = avma;
606 7 : T = RgX_to_FpX(T,p);
607 7 : x = RgX_to_FpX(x, p);
608 7 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
609 7 : z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
610 7 : return gc_upto(av, z);
611 : }
612 182 : case t_INT:
613 182 : case t_FRAC: break;
614 0 : default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
615 : }
616 : }
617 49 : if (typ(x) != t_POL)
618 : {
619 0 : if (!is_rational_t(typ(x)))
620 0 : pari_err_TYPE("conjvec [not a rational t_POL]",x);
621 0 : retconst_col(lx-3, gcopy(x));
622 : }
623 49 : RgX_check_QX(x,"conjvec");
624 49 : av = avma;
625 49 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
626 49 : r = cleanroots(T,prec);
627 49 : z = cgetg(lx-2,t_COL);
628 182 : for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
629 49 : return gc_upto(av, z);
630 : }
631 :
632 0 : default:
633 0 : pari_err_TYPE("conjvec",x);
634 : return NULL; /* LCOV_EXCL_LINE */
635 : }
636 0 : return z;
637 : }
638 :
639 : /********************************************************************/
640 : /** **/
641 : /** ADDITION **/
642 : /** **/
643 : /********************************************************************/
644 : static GEN
645 24594886 : mkpadic_mod(GEN u, GEN p, GEN pd, long e, long d)
646 24594886 : { retmkpadic(modii(u, pd), icopy(p), icopy(pd), e, d); }
647 :
648 : /* x, y compatible PADIC, op = add or sub */
649 : static GEN
650 19974791 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
651 : {
652 19974791 : pari_sp av = avma;
653 : long d, e, r, rx, ry;
654 19974791 : GEN u, z, mod, pdx, pdy, ux, uy, p = padic_p(x);
655 : int swap;
656 :
657 19974791 : e = valp(x);
658 19974791 : r = valp(y); d = r-e;
659 19974791 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
660 19974791 : pdx = padic_pd(x); ux = padic_u(x);
661 19974791 : pdy = padic_pd(y); uy = padic_u(y);
662 19974791 : rx = precp(x);
663 19974791 : ry = precp(y);
664 19974791 : if (d) /* v(x) < v(y) */
665 : {
666 10704514 : r = d+ry; z = powiu(p,d);
667 10704621 : if (r < rx) mod = mulii(z,pdy); else { r = rx; mod = pdx; }
668 10704639 : z = mulii(z, uy);
669 10704588 : u = swap? op(z, ux): op(ux, z);
670 : }
671 : else
672 : {
673 : long c;
674 9270277 : if (ry < rx) { r=ry; mod = pdy; } else { r=rx; mod = pdx; }
675 9270277 : u = swap? op(uy, ux): op(ux, uy);
676 9270910 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
677 : {
678 138875 : set_avma(av); return zeropadic(p, e+r);
679 : }
680 9131874 : if (c)
681 : {
682 3351658 : mod = diviiexact(mod, powiu(p,c));
683 3351660 : r -= c;
684 3351660 : e += c;
685 : }
686 : }
687 19836436 : return gc_upto(av, mkpadic_mod(u, p, mod, e, r));
688 : }
689 : /* Rg_to_Fp(t_FRAC) without GC */
690 : static GEN
691 226838 : Q_to_Fp(GEN x, GEN p)
692 226838 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
693 : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
694 : static GEN
695 4762976 : addQp(GEN x, GEN y)
696 : {
697 4762976 : pari_sp av = avma;
698 4762976 : long d, r, e, vy = valp(y), py = precp(y);
699 4762976 : GEN q, mod, u, p = padic_p(y);
700 :
701 4762976 : e = Q_pvalrem(x, p, &x);
702 4762974 : d = vy - e; r = d + py;
703 4762974 : if (r <= 0) { set_avma(av); return gcopy(y); }
704 4761111 : mod = padic_pd(y);
705 4761111 : u = padic_u(y);
706 4761111 : if (d > 0)
707 : {
708 1376568 : q = powiu(p,d);
709 1376572 : mod = mulii(mod, q);
710 1376569 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
711 1376569 : u = addii(x, mulii(u, q));
712 : }
713 3384543 : else if (d < 0)
714 : {
715 405332 : q = powiu(p,-d);
716 405332 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
717 405332 : u = addii(u, mulii(x, q));
718 405332 : r = py; e = vy;
719 : }
720 : else
721 : {
722 : long c;
723 2979211 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
724 2979211 : u = addii(u, x);
725 2979211 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
726 : {
727 1204 : set_avma(av); return zeropadic(p,e+r);
728 : }
729 2977992 : if (c)
730 : {
731 970686 : mod = diviiexact(mod, powiu(p,c));
732 970686 : r -= c;
733 970686 : e += c;
734 : }
735 : }
736 4759893 : return gc_upto(av, mkpadic_mod(u, p, mod, e, r));
737 : }
738 :
739 : /* Mod(x,X) + Mod(y,X) */
740 : #define addsub_polmod_same addsub_polmod_scal
741 : /* Mod(x,X) +/- Mod(y,Y) */
742 : static GEN
743 7203 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
744 : {
745 7203 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
746 7203 : GEN z = cgetg(3,t_POLMOD);
747 7203 : long vx = varn(X), vy = varn(Y);
748 7203 : if (vx==vy) {
749 : pari_sp av;
750 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
751 14 : warn_coercion(X,Y,gel(z,1));
752 14 : gel(z,2) = gc_upto(av, gmod(op(x, y), gel(z,1))); return z;
753 : }
754 7189 : if (varncmp(vx, vy) < 0)
755 7189 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
756 : else
757 0 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
758 7189 : gel(z,2) = op(x, y); return z;
759 : }
760 : /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
761 : static GEN
762 13456529 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
763 13456529 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
764 :
765 : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
766 : static GEN
767 414436 : add_ser_scal(GEN y, GEN x)
768 : {
769 : long i, v, ly, vy;
770 : GEN z;
771 :
772 414436 : if (isrationalzero(x)) return gcopy(y);
773 387675 : ly = lg(y);
774 387675 : v = valser(y);
775 387675 : if (v < 3-ly) return gcopy(y);
776 : /* v + ly >= 3 */
777 387395 : if (v < 0)
778 : {
779 1176 : z = cgetg(ly,t_SER); z[1] = y[1];
780 3304 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
781 1176 : gel(z,i) = gadd(x,gel(y,i)); i++;
782 3087 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
783 1176 : return normalizeser(z);
784 : }
785 386219 : vy = varn(y);
786 386219 : if (v > 0)
787 : {
788 19649 : if (ser_isexactzero(y))
789 9485 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
790 10164 : y -= v; ly += v;
791 10164 : z = cgetg(ly,t_SER);
792 10164 : x = gcopy(x);
793 20762 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
794 : }
795 366570 : else if (ser_isexactzero(y)) return gcopy(y);
796 : else
797 : { /* v = 0, ly >= 3 */
798 366563 : z = cgetg(ly,t_SER);
799 366563 : x = gadd(x, gel(y,2));
800 366563 : i = 3;
801 : }
802 1613567 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
803 376727 : gel(z,2) = x;
804 376727 : z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
805 376727 : return gequal0(x)? normalizeser(z): z;
806 : }
807 : static long
808 7185086 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
809 : /* x,y t_SER in the same variable: x+y */
810 : static GEN
811 3592942 : ser_add(GEN x, GEN y)
812 : {
813 3592942 : long i, lx,ly, n = valser(y) - valser(x);
814 : GEN z;
815 3592942 : if (n < 0) { n = -n; swap(x,y); }
816 : /* valser(x) <= valser(y) */
817 3592942 : lx = _serprec(x);
818 3592942 : if (lx == 2) /* don't lose type information */
819 : {
820 798 : z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
821 798 : setvalser(z, valser(x)); return z;
822 : }
823 3592144 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
824 3592144 : if (n)
825 : {
826 107831 : if (n+2 > lx) return gcopy(x);
827 107131 : z = cgetg(ly,t_SER);
828 819284 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
829 506743 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
830 : } else {
831 3484313 : z = cgetg(ly,t_SER);
832 17642238 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
833 : }
834 3591444 : z[1] = x[1]; return normalizeser(z);
835 : }
836 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
837 : static GEN
838 8812213 : add_rfrac_scal(GEN y, GEN x)
839 : {
840 : pari_sp av;
841 : GEN n;
842 :
843 8812213 : if (isintzero(x)) return gcopy(y); /* frequent special case */
844 5138455 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
845 5138455 : return gc_upto(av, gred_rfrac_simple(n, gel(y,2)));
846 : }
847 :
848 : /* x "scalar", ty != t_MAT and nonscalar */
849 : static GEN
850 42826811 : add_scal(GEN y, GEN x, long ty)
851 : {
852 42826811 : switch(ty)
853 : {
854 37886542 : case t_POL: return RgX_Rg_add(y, x);
855 414408 : case t_SER: return add_ser_scal(y, x);
856 4480454 : case t_RFRAC: return add_rfrac_scal(y, x);
857 0 : case t_COL: return RgC_Rg_add(y, x);
858 45395 : case t_VEC:
859 45395 : if (isintzero(x)) return gcopy(y);
860 175 : break;
861 : }
862 187 : pari_err_TYPE2("+",x,y);
863 : return NULL; /* LCOV_EXCL_LINE */
864 : }
865 :
866 : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
867 : static GEN
868 14999463 : setfrac(GEN z, GEN a, GEN b)
869 : {
870 14999463 : gel(z,1) = icopy_avma(a, (pari_sp)z);
871 14999463 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
872 14999523 : set_avma((pari_sp)gel(z,2)); return z;
873 : }
874 : /* z <- a / (b*Q), (Q,a) = 1 */
875 : static GEN
876 14120656 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
877 : {
878 14120656 : GEN q = Qdivii(a, b); /* != 0 */
879 14120723 : if (typ(q) == t_INT)
880 : {
881 1904476 : gel(z,1) = gc_INT((pari_sp)Q, q);
882 1904476 : gel(z,2) = Q; return z;
883 : }
884 12216247 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
885 : }
886 : static GEN
887 40598719 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
888 : {
889 40598719 : GEN x1 = gel(x,1), x2 = gel(x,2);
890 40598719 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
891 40598719 : int s = cmpii(x2, y2);
892 :
893 : /* common denominator: (x1 op y1) / x2 */
894 40598711 : if (!s)
895 : {
896 15299033 : pari_sp av = avma;
897 15299033 : return gc_upto(av, Qdivii(op(x1, y1), x2));
898 : }
899 25299678 : z = cgetg(3, t_FRAC);
900 25299696 : if (s < 0)
901 : {
902 17684778 : Q = dvmdii(y2, x2, &r);
903 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
904 17684742 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
905 9722624 : delta = gcdii(x2,r);
906 : }
907 : else
908 : {
909 7614918 : Q = dvmdii(x2, y2, &r);
910 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
911 7614935 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
912 1456328 : delta = gcdii(y2,r);
913 : }
914 : /* delta = gcd(x2,y2) */
915 11178978 : if (equali1(delta))
916 : { /* numerator is nonzero */
917 8395698 : gel(z,1) = gc_INT((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
918 8395700 : gel(z,2) = mulii(x2,y2); return z;
919 : }
920 2783275 : x2 = diviiexact(x2,delta);
921 2783280 : y2 = diviiexact(y2,delta);
922 2783279 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
923 2783281 : q = dvmdii(n, delta, &r);
924 2783279 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
925 2529254 : r = gcdii(delta, r);
926 2529258 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
927 2529258 : return setfrac(z, n, mulii(mulii(x2, y2), delta));
928 : }
929 :
930 : /* assume x2, y2 are t_POLs in the same variable */
931 : static GEN
932 3017228 : add_rfrac(GEN x, GEN y)
933 : {
934 3017228 : pari_sp av = avma;
935 3017228 : GEN x1 = gel(x,1), x2 = gel(x,2);
936 3017228 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
937 :
938 3017228 : delta = RgX_gcd(x2,y2);
939 3017228 : if (!degpol(delta))
940 : {
941 644 : n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
942 644 : d = RgX_mul(x2, y2);
943 644 : return gc_upto(av, gred_rfrac_simple(n, d));
944 : }
945 3016584 : x2 = RgX_div(x2,delta);
946 3016584 : y2 = RgX_div(y2,delta);
947 3016584 : n = gadd(gmul(x1,y2), gmul(y1,x2));
948 3016584 : if (!signe(n))
949 : {
950 721418 : n = simplify_shallow(n);
951 721418 : if (isexactzero(n))
952 : {
953 721411 : if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
954 7 : return gc_upto(av, scalarpol(n, varn(x2)));
955 : }
956 7 : return gc_GEN(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
957 : }
958 2295166 : if (degpol(n) == 0)
959 1150589 : return gc_upto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
960 1144577 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
961 1144577 : if (isexactzero(r))
962 : {
963 : GEN z;
964 227996 : d = RgX_mul(x2, y2);
965 : /* "constant" denominator ? */
966 227996 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
967 227996 : return gc_upto(av, z);
968 : }
969 916581 : r = RgX_gcd(delta, r);
970 916581 : if (degpol(r))
971 : {
972 160593 : n = RgX_div(n, r);
973 160593 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
974 : }
975 : else
976 755988 : d = RgX_mul(gel(x,2), y2);
977 916581 : return gc_upto(av, gred_rfrac_simple(n, d));
978 : }
979 :
980 : GEN
981 6149360256 : gadd(GEN x, GEN y)
982 : {
983 6149360256 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
984 : pari_sp av;
985 : GEN z, p1;
986 :
987 6149360256 : if (tx == ty) switch(tx) /* shortcut to generic case */
988 : {
989 3061886080 : case t_INT: return addii(x,y);
990 1659439615 : case t_REAL: return addrr(x,y);
991 1170566 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
992 1170566 : z = cgetg(3,t_INTMOD);
993 1170566 : if (X==Y || equalii(X,Y))
994 1170552 : return add_intmod_same(z, X, gel(x,2), gel(y,2));
995 14 : gel(z,1) = gcdii(X,Y);
996 14 : warn_coercion(X,Y,gel(z,1));
997 14 : av = avma; p1 = addii(gel(x,2),gel(y,2));
998 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
999 : }
1000 34505332 : case t_FRAC: return addsub_frac(x,y,addii);
1001 360806977 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1002 360335994 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1003 360853381 : if (isintzero(gel(z,2)))
1004 : {
1005 521269 : set_avma((pari_sp)(z+3));
1006 521269 : return gadd(gel(x,1),gel(y,1));
1007 : }
1008 360193263 : gel(z,1) = gadd(gel(x,1),gel(y,1));
1009 360286392 : return z;
1010 14532634 : case t_PADIC:
1011 14532634 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1012 14532259 : return addsub_pp(x,y, addii);
1013 672 : case t_QUAD: z = cgetg(4,t_QUAD);
1014 672 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1015 672 : gel(z,1) = ZX_copy(gel(x,1));
1016 672 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1017 672 : gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
1018 10151602 : case t_POLMOD:
1019 10151602 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1020 10144434 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
1021 7168 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
1022 13133901 : case t_FFELT: return FF_add(x,y);
1023 138939362 : case t_POL:
1024 138939362 : vx = varn(x);
1025 138939362 : vy = varn(y);
1026 138939362 : if (vx != vy) {
1027 4863983 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1028 1801550 : else return RgX_Rg_add(y, x);
1029 : }
1030 134075379 : return RgX_add(x, y);
1031 3587580 : case t_SER:
1032 3587580 : vx = varn(x);
1033 3587580 : vy = varn(y);
1034 3587580 : if (vx != vy) {
1035 28 : if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1036 21 : else return add_ser_scal(y, x);
1037 : }
1038 3587552 : return ser_add(x, y);
1039 4309129 : case t_RFRAC:
1040 4309129 : vx = varn(gel(x,2));
1041 4309129 : vy = varn(gel(y,2));
1042 4309129 : if (vx != vy) {
1043 1291901 : if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1044 538397 : else return add_rfrac_scal(y, x);
1045 : }
1046 3017228 : return add_rfrac(x,y);
1047 2115710 : case t_VEC:
1048 2115710 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1049 2115710 : return RgV_add(x,y);
1050 1113425 : case t_COL:
1051 1113425 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1052 1113425 : return RgC_add(x,y);
1053 671412 : case t_MAT:
1054 671412 : lx = lg(x);
1055 671412 : if (lg(y) != lx) pari_err_OP("+",x,y);
1056 671412 : if (lx == 1) return cgetg(1, t_MAT);
1057 671412 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1058 671405 : return RgM_add(x,y);
1059 :
1060 0 : default: pari_err_TYPE2("+",x,y);
1061 : }
1062 : /* tx != ty */
1063 846714198 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1064 :
1065 846714198 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1066 : {
1067 732826506 : case t_INT:
1068 : switch(ty)
1069 : {
1070 442826807 : case t_REAL: return addir(x,y);
1071 2138717 : case t_INTMOD:
1072 2138717 : z = cgetg(3, t_INTMOD);
1073 2138717 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1074 27566402 : case t_FRAC: z = cgetg(3,t_FRAC);
1075 27566404 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1076 27566365 : gel(z,2) = icopy(gel(y,2)); return z;
1077 253371843 : case t_COMPLEX: return addRc(x, y);
1078 5596984 : case t_PADIC:
1079 5596984 : if (!signe(x)) return gcopy(y);
1080 4533203 : return addQp(x,y);
1081 1113 : case t_QUAD: return addRq(x, y);
1082 1365028 : case t_FFELT: return FF_Z_add(y,x);
1083 : }
1084 :
1085 : case t_REAL:
1086 58623846 : switch(ty)
1087 : {
1088 14406373 : case t_FRAC:
1089 14406373 : if (!signe(gel(y,1))) return rcopy(x);
1090 14406373 : if (!signe(x))
1091 : {
1092 12080 : lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1093 12080 : return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1094 : }
1095 14394293 : av = avma; z = addir(gel(y,1), mulir(gel(y,2),x));
1096 14392915 : return gc_leaf(av, divri(z,gel(y,2)));
1097 44217403 : case t_COMPLEX: return addRc(x, y);
1098 70 : case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, realprec(x));
1099 :
1100 0 : default: pari_err_TYPE2("+",x,y);
1101 : }
1102 :
1103 17647 : case t_INTMOD:
1104 : switch(ty)
1105 : {
1106 17507 : case t_FRAC: { GEN X = gel(x,1);
1107 17507 : z = cgetg(3, t_INTMOD);
1108 17507 : p1 = Fp_div(gel(y,1), gel(y,2), X);
1109 17507 : return add_intmod_same(z, X, p1, gel(x,2));
1110 : }
1111 14 : case t_FFELT:
1112 14 : if (!equalii(gel(x,1),FF_p_i(y)))
1113 0 : pari_err_OP("+",x,y);
1114 14 : return FF_Z_add(y,gel(x,2));
1115 91 : case t_COMPLEX: return addRc(x, y);
1116 0 : case t_PADIC: { GEN X = gel(x,1);
1117 0 : z = cgetg(3, t_INTMOD);
1118 0 : return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1119 : }
1120 35 : case t_QUAD: return addRq(x, y);
1121 : }
1122 :
1123 : case t_FRAC:
1124 : switch (ty)
1125 : {
1126 10327658 : case t_COMPLEX: return addRc(x, y);
1127 229771 : case t_PADIC:
1128 229771 : if (!signe(gel(x,1))) return gcopy(y);
1129 229771 : return addQp(x,y);
1130 133 : case t_QUAD: return addRq(x, y);
1131 1337 : case t_FFELT: return FF_Q_add(y, x);
1132 : }
1133 :
1134 : case t_FFELT:
1135 1 : pari_err_TYPE2("+",x,y);
1136 :
1137 35 : case t_COMPLEX:
1138 : switch(ty)
1139 : {
1140 28 : case t_PADIC:
1141 28 : return Zp_nosquare_m1(padic_p(y))? addRc(y, x): addTp(x, y);
1142 7 : case t_QUAD:
1143 7 : lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1144 7 : return gequal0(y)? gcopy(x): addqf(y, x, lx);
1145 : }
1146 :
1147 : case t_PADIC: /* ty == t_QUAD */
1148 6 : return (kro_quad(y, padic_p(x)) == -1)? addRq(x, y): addTp(y, x);
1149 : }
1150 : /* tx < ty, !is_const_t(y) */
1151 48507075 : switch(ty)
1152 : {
1153 7688 : case t_MAT:
1154 7688 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1155 7688 : if (isrationalzero(x)) return gcopy(y);
1156 7611 : return RgM_Rg_add(y, x);
1157 186359 : case t_COL:
1158 186359 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1159 186359 : return RgC_Rg_add(y, x);
1160 2409837 : case t_POLMOD: /* is_const_t(tx) in this case */
1161 2409837 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1162 : }
1163 45903191 : if (is_scalar_t(tx)) {
1164 42834654 : if (tx == t_POLMOD)
1165 : {
1166 117831 : vx = varn(gel(x,1));
1167 117831 : vy = gvar(y);
1168 117831 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1169 : else
1170 87780 : if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1171 30590 : return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1172 : }
1173 42716823 : return add_scal(y, x, ty);
1174 : }
1175 : /* x and y are not scalars, ty != t_MAT */
1176 3068530 : vx = gvar(x);
1177 3068532 : vy = gvar(y);
1178 3068532 : if (vx != vy) { /* x or y is treated as a scalar */
1179 22759 : if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1180 32417 : return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1181 32417 : : add_scal(y, x, ty);
1182 : }
1183 : /* vx = vy */
1184 3045773 : switch(tx)
1185 : {
1186 3045276 : case t_POL:
1187 : switch (ty)
1188 : {
1189 5418 : case t_SER:
1190 5418 : if (lg(x) == 2) return gcopy(y);
1191 5397 : i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1192 5397 : i = lg(y) + valser(y) - i;
1193 5397 : if (i < 3) return gcopy(y);
1194 5390 : p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1195 5390 : settyp(p1, t_VECSMALL); /* p1 left on stack */
1196 5390 : return y;
1197 :
1198 3039858 : case t_RFRAC: return add_rfrac_scal(y, x);
1199 : }
1200 0 : break;
1201 :
1202 497 : case t_SER:
1203 497 : if (ty == t_RFRAC)
1204 : {
1205 : long vn, vd;
1206 497 : av = avma;
1207 497 : vn = gval(gel(y,1), vy);
1208 497 : vd = RgX_valrem_inexact(gel(y,2), NULL);
1209 497 : if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
1210 :
1211 497 : l = lg(x) + valser(x) - (vn - vd);
1212 497 : if (l < 3) { set_avma(av); return gcopy(x); }
1213 497 : return gc_upto(av, gadd(x, rfrac_to_ser_i(y, l)));
1214 : }
1215 0 : break;
1216 : }
1217 0 : pari_err_TYPE2("+",x,y);
1218 : return NULL; /* LCOV_EXCL_LINE */
1219 : }
1220 :
1221 : GEN
1222 270740483 : gaddsg(long x, GEN y)
1223 : {
1224 270740483 : long ty = typ(y);
1225 : GEN z;
1226 :
1227 270740483 : switch(ty)
1228 : {
1229 124845352 : case t_INT: return addsi(x,y);
1230 119267190 : case t_REAL: return addsr(x,y);
1231 12059 : case t_INTMOD:
1232 12059 : z = cgetg(3, t_INTMOD);
1233 12059 : return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1234 15250631 : case t_FRAC: z = cgetg(3,t_FRAC);
1235 15250631 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1236 15250631 : gel(z,2) = icopy(gel(y,2)); return z;
1237 8226730 : case t_COMPLEX:
1238 8226730 : z = cgetg(3, t_COMPLEX);
1239 8226730 : gel(z,1) = gaddsg(x, gel(y,1));
1240 8226730 : gel(z,2) = gcopy(gel(y,2)); return z;
1241 :
1242 3138521 : default: return gadd(stoi(x), y);
1243 : }
1244 : }
1245 :
1246 : GEN
1247 3235140 : gsubsg(long x, GEN y)
1248 : {
1249 : GEN z, a, b;
1250 : pari_sp av;
1251 :
1252 3235140 : switch(typ(y))
1253 : {
1254 276974 : case t_INT: return subsi(x,y);
1255 1230556 : case t_REAL: return subsr(x,y);
1256 56 : case t_INTMOD:
1257 56 : z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1258 56 : return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1259 732542 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1260 732542 : gel(z,1) = gc_INT((pari_sp)z, subii(mulis(b,x), a));
1261 732542 : gel(z,2) = icopy(gel(y,2)); return z;
1262 957856 : case t_COMPLEX:
1263 957856 : z = cgetg(3, t_COMPLEX);
1264 957856 : gel(z,1) = gsubsg(x, gel(y,1));
1265 957856 : gel(z,2) = gneg(gel(y,2)); return z;
1266 : }
1267 37156 : av = avma;
1268 37156 : return gc_upto(av, gadd(stoi(x), gneg_i(y)));
1269 : }
1270 :
1271 : /********************************************************************/
1272 : /** **/
1273 : /** SUBTRACTION **/
1274 : /** **/
1275 : /********************************************************************/
1276 :
1277 : GEN
1278 2887715683 : gsub(GEN x, GEN y)
1279 : {
1280 2887715683 : long tx = typ(x), ty = typ(y);
1281 : pari_sp av;
1282 : GEN z;
1283 2887715683 : if (tx == ty) switch(tx) /* shortcut to generic case */
1284 : {
1285 2094523118 : case t_INT: return subii(x,y);
1286 586890880 : case t_REAL: return subrr(x,y);
1287 1021783 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1288 1021783 : z = cgetg(3,t_INTMOD);
1289 1021783 : if (X==Y || equalii(X,Y))
1290 1021769 : return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1291 14 : gel(z,1) = gcdii(X,Y);
1292 14 : warn_coercion(X,Y,gel(z,1));
1293 14 : av = avma; p1 = subii(gel(x,2),gel(y,2));
1294 14 : gel(z,2) = gc_INT(av, modii(p1, gel(z,1))); return z;
1295 : }
1296 6093398 : case t_FRAC: return addsub_frac(x,y, subii);
1297 103069626 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1298 103100708 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1299 103021319 : if (isintzero(gel(z,2)))
1300 : {
1301 21322 : set_avma((pari_sp)(z+3));
1302 21322 : return gsub(gel(x,1),gel(y,1));
1303 : }
1304 103003482 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1305 102990645 : return z;
1306 5442581 : case t_PADIC:
1307 5442581 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1308 5442581 : return addsub_pp(x,y, subii);
1309 168 : case t_QUAD: z = cgetg(4,t_QUAD);
1310 168 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1311 168 : gel(z,1) = ZX_copy(gel(x,1));
1312 168 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1313 168 : gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1314 871703 : case t_POLMOD:
1315 871703 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1316 871668 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1317 35 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1318 240344 : case t_FFELT: return FF_sub(x,y);
1319 20216032 : case t_POL: {
1320 20216032 : long vx = varn(x);
1321 20216032 : long vy = varn(y);
1322 20216032 : if (vx != vy) {
1323 1150623 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1324 85352 : else return Rg_RgX_sub(x, y);
1325 : }
1326 19065409 : return RgX_sub(x, y);
1327 : }
1328 299580 : case t_VEC:
1329 299580 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1330 299580 : return RgV_sub(x,y);
1331 3604746 : case t_COL:
1332 3604746 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1333 3604746 : return RgC_sub(x,y);
1334 254334 : case t_MAT: {
1335 254334 : long lx = lg(x);
1336 254334 : if (lg(y) != lx) pari_err_OP("+",x,y);
1337 254337 : if (lx == 1) return cgetg(1, t_MAT);
1338 174231 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1339 174231 : return RgM_sub(x,y);
1340 : }
1341 2449645 : case t_RFRAC: case t_SER: break;
1342 :
1343 0 : default: pari_err_TYPE2("+",x,y);
1344 : }
1345 66581433 : av = avma;
1346 66581433 : return gc_upto(av, gadd(x,gneg_i(y)));
1347 : }
1348 :
1349 : /********************************************************************/
1350 : /** **/
1351 : /** MULTIPLICATION **/
1352 : /** **/
1353 : /********************************************************************/
1354 : static GEN
1355 321174 : mul_ser_scal(GEN x, GEN t)
1356 : {
1357 321174 : if (isexactzero(t)) return gmul(Rg_get_0(x), t);
1358 317926 : if (isint1(t)) return gcopy(x);
1359 265615 : if (ser_isexactzero(x))
1360 : {
1361 378 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gmul(gel(x,2), t), varn(x), 1);
1362 378 : setvalser(z, valser(x)); return z;
1363 : }
1364 3454570 : pari_APPLY_ser(gmul(gel(x,i), t));
1365 : }
1366 : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1367 : * [n/d a valid RFRAC] */
1368 : static GEN
1369 10454093 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1370 : {
1371 10454093 : pari_sp av = avma;
1372 : GEN z;
1373 :
1374 10454093 : switch(typ(x))
1375 : {
1376 21 : case t_PADIC:
1377 21 : n = gmul(n, x);
1378 21 : d = gcvtop(d, padic_p(x), signe(padic_u(x))? precp(x): 1);
1379 21 : return gc_upto(av, gdiv(n,d));
1380 :
1381 938 : case t_INTMOD: case t_POLMOD:
1382 938 : n = gmul(n, x);
1383 938 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1384 938 : return gc_upto(av, gdiv(n,d));
1385 : }
1386 10453134 : z = gred_rfrac2(x, d);
1387 10453134 : n = simplify_shallow(n);
1388 10453134 : if (typ(z) == t_RFRAC)
1389 : {
1390 7928497 : n = gmul(gel(z,1), n);
1391 7928497 : d = gel(z,2);
1392 7928497 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1393 0 : z = RgX_Rg_div(n, d);
1394 : else
1395 7928497 : z = gred_rfrac_simple(n, d);
1396 : }
1397 : else
1398 2524637 : z = gmul(z, n);
1399 10453134 : return gc_upto(av, z);
1400 : }
1401 : static GEN
1402 154346028 : mul_scal(GEN y, GEN x, long ty)
1403 : {
1404 154346028 : switch(ty)
1405 : {
1406 145617804 : case t_POL:
1407 145617804 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1408 133552353 : return RgX_Rg_mul(y, x);
1409 313082 : case t_SER: return mul_ser_scal(y, x);
1410 8415122 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1411 14 : case t_QFB:
1412 14 : if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1413 : }
1414 12 : pari_err_TYPE2("*",x,y);
1415 : return NULL; /* LCOV_EXCL_LINE */
1416 : }
1417 : static GEN
1418 7645958 : mul_self_scal(GEN x, GEN y)
1419 641090731 : { pari_APPLY_same(gmul(y,gel(x,i))); }
1420 :
1421 : static GEN
1422 161432 : mul_gen_rfrac(GEN X, GEN Y)
1423 : {
1424 161432 : GEN y1 = gel(Y,1), y2 = gel(Y,2);
1425 161432 : long vx = gvar(X), vy = varn(y2);
1426 167193 : return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1427 5761 : gred_rfrac_simple(gmul(y1,X), y2);
1428 : }
1429 : /* (x1/x2) * (y1/y2) */
1430 : static GEN
1431 7910299 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1432 : {
1433 : GEN z, X, Y;
1434 7910299 : pari_sp av = avma;
1435 :
1436 7910299 : X = gred_rfrac2(x1, y2);
1437 7910299 : Y = gred_rfrac2(y1, x2);
1438 7910299 : if (typ(X) == t_RFRAC)
1439 : {
1440 6631203 : if (typ(Y) == t_RFRAC) {
1441 6565046 : x1 = gel(X,1);
1442 6565046 : x2 = gel(X,2);
1443 6565046 : y1 = gel(Y,1);
1444 6565046 : y2 = gel(Y,2);
1445 6565046 : z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1446 : } else
1447 66157 : z = mul_gen_rfrac(Y, X);
1448 : }
1449 1279096 : else if (typ(Y) == t_RFRAC)
1450 95275 : z = mul_gen_rfrac(X, Y);
1451 : else
1452 1183821 : z = gmul(X, Y);
1453 7910299 : return gc_upto(av, z);
1454 : }
1455 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1456 : static GEN
1457 262365 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1458 : {
1459 262365 : pari_sp av = avma;
1460 262365 : GEN X = gred_rfrac2(x1, y2);
1461 262365 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1462 : {
1463 255785 : x2 = RgX_mul(gel(X,2), x2);
1464 255785 : x1 = gel(X,1);
1465 : }
1466 : else
1467 6580 : x1 = X;
1468 262365 : return gc_upto(av, gred_rfrac_simple(x1, x2));
1469 : }
1470 :
1471 : /* Mod(y, Y) * x, assuming x scalar */
1472 : static GEN
1473 3238992 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1474 3238992 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1475 :
1476 : /* cf mulqq */
1477 : static GEN
1478 5860209 : quad_polmod_mul(GEN T, GEN x, GEN y)
1479 : {
1480 5860209 : GEN b = gel(T,3), c = gel(T,2);
1481 5860209 : GEN ux = gel(x,2), vx = gel(x,3), uy = gel(y,2), vy = gel(y,3);
1482 : GEN z, s, U, V, E, F;
1483 : pari_sp av, av2;
1484 5860209 : z = cgetg(4, t_POL); z[1] = x[1]; av = avma;
1485 5860209 : U = gmul(ux, uy);
1486 5860209 : V = gmul(vx, vy); s = gmul(c, V);
1487 5860209 : E = gmul(gadd(ux, vx), gadd(uy, vy));
1488 5860209 : if (typ(b) != t_INT) F = gadd(U, gmul(gaddgs(b, 1), V)); /* = U + (b+1) V */
1489 : else
1490 : { /* minor optimization */
1491 5860188 : if (!signe(b)) F = gadd(U, V); /* b = 0 */
1492 4284989 : else if (is_pm1(b)) /* b = 1 or -1 */
1493 4284331 : F = signe(b) < 0? U: gadd(U, gmul2n(V,1));
1494 : else /* |b| > 1 */
1495 658 : F = gadd(U, gmul(addis(b, 1), V));
1496 : }
1497 5860209 : av2 = avma;
1498 5860209 : gel(z,2) = gsub(U, s);
1499 5860209 : gel(z,3) = gsub(E, F);
1500 5860209 : gc_slice_unsafe(av, av2, z+2, 2);
1501 5860209 : return normalizepol_lg(z,4);
1502 : }
1503 : /* Mod(x,T) * Mod(y,T) */
1504 : static GEN
1505 8529774 : mul_polmod_same(GEN T, GEN x, GEN y)
1506 : {
1507 8529774 : GEN z = cgetg(3,t_POLMOD), a;
1508 8529774 : long v = varn(T), lx = lg(x), ly = lg(y);
1509 8529774 : gel(z,1) = RgX_copy(T);
1510 : /* x * y mod T optimised */
1511 8529774 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1512 7902201 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1513 1546705 : a = gmul(x, y);
1514 : else
1515 : {
1516 6983069 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1517 5854931 : a = quad_polmod_mul(T, x, y);
1518 : else
1519 1128138 : a = RgXQ_mul(x, y, gel(z,1));
1520 : }
1521 8529774 : gel(z,2) = a; return z;
1522 : }
1523 : static GEN
1524 419385 : sqr_polmod(GEN T, GEN x)
1525 : {
1526 419385 : GEN a, z = cgetg(3,t_POLMOD);
1527 419385 : gel(z,1) = RgX_copy(T);
1528 419385 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1529 28967 : a = gsqr(x);
1530 : else
1531 : {
1532 390418 : pari_sp av = avma;
1533 390418 : a = RgXQ_sqr(x, gel(z,1));
1534 390418 : a = gc_upto(av, a);
1535 : }
1536 419385 : gel(z,2) = a; return z;
1537 : }
1538 : /* Mod(x,X) * Mod(y,Y) */
1539 : static GEN
1540 2674945 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1541 : {
1542 2674945 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1543 2674945 : long vx = varn(X), vy = varn(Y);
1544 2674945 : GEN z = cgetg(3,t_POLMOD);
1545 :
1546 2674945 : if (vx==vy) {
1547 : pari_sp av;
1548 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
1549 14 : warn_coercion(X,Y,gel(z,1));
1550 14 : gel(z,2) = gc_upto(av, gmod(gmul(x, y), gel(z,1)));
1551 14 : return z;
1552 : }
1553 2674931 : if (varncmp(vx, vy) < 0)
1554 410662 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1555 : else
1556 2264269 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1557 2674931 : gel(z,2) = gmul(x, y); return z;
1558 : }
1559 :
1560 : #if 0 /* used by 3M only */
1561 : /* set z = x+y and return 1 if x,y have the same sign
1562 : * set z = x-y and return 0 otherwise */
1563 : static int
1564 : did_add(GEN x, GEN y, GEN *z)
1565 : {
1566 : long tx = typ(x), ty = typ(y);
1567 : if (tx == ty) switch(tx)
1568 : {
1569 : case t_INT: *z = addii(x,y); return 1;
1570 : case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1571 : case t_REAL:
1572 : if (signe(x) == -signe(y))
1573 : { *z = subrr(x,y); return 0; }
1574 : else
1575 : { *z = addrr(x,y); return 1; }
1576 : }
1577 : if (tx == t_REAL) switch(ty)
1578 : {
1579 : case t_INT:
1580 : if (signe(x) == -signe(y))
1581 : { *z = subri(x,y); return 0; }
1582 : else
1583 : { *z = addri(x,y); return 1; }
1584 : case t_FRAC:
1585 : if (signe(x) == -signe(gel(y,1)))
1586 : { *z = gsub(x,y); return 0; }
1587 : else
1588 : { *z = gadd(x,y); return 1; }
1589 : }
1590 : else if (ty == t_REAL) switch(tx)
1591 : {
1592 : case t_INT:
1593 : if (signe(x) == -signe(y))
1594 : { *z = subir(x,y); return 0; }
1595 : else
1596 : { *z = addir(x,y); return 1; }
1597 : case t_FRAC:
1598 : if (signe(gel(x,1)) == -signe(y))
1599 : { *z = gsub(x,y); return 0; }
1600 : else
1601 : { *z = gadd(x,y); return 1; }
1602 : }
1603 : *z = gadd(x,y); return 1;
1604 : }
1605 : #endif
1606 : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1607 : static GEN
1608 11646834 : mulcIR(GEN x, GEN y)
1609 : {
1610 11646834 : GEN z = cgetg(3,t_COMPLEX);
1611 11646992 : pari_sp av = avma;
1612 11646992 : gel(z,1) = gc_upto(av, gneg(gmul(y,gel(x,2))));
1613 11647134 : gel(z,2) = gmul(y, gel(x,1));
1614 11647112 : return z;
1615 :
1616 : }
1617 : /* x,y COMPLEX */
1618 : static GEN
1619 279036320 : mulcc(GEN x, GEN y)
1620 : {
1621 279036320 : GEN xr = gel(x,1), xi = gel(x,2);
1622 279036320 : GEN yr = gel(y,1), yi = gel(y,2);
1623 : GEN p1, p2, p3, p4, z;
1624 : pari_sp tetpil, av;
1625 :
1626 279036320 : if (isintzero(xr))
1627 : {
1628 15644925 : if (isintzero(yr)) {
1629 7122203 : av = avma;
1630 7122203 : return gc_upto(av, gneg(gmul(xi,yi)));
1631 : }
1632 8522571 : return mulcIR(y, xi);
1633 : }
1634 263385075 : if (isintzero(yr)) return mulcIR(x, yi);
1635 :
1636 260265726 : z = cgetg(3,t_COMPLEX); av = avma;
1637 : #if 0
1638 : /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1639 : * e.g. xr + xi if exponents differ */
1640 : if (did_add(xr, xi, &p3))
1641 : {
1642 : if (did_add(yr, yi, &p4)) {
1643 : /* R = xr*yr - xi*yi
1644 : * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1645 : p1 = gmul(xr,yr);
1646 : p2 = gmul(xi,yi); p2 = gneg(p2);
1647 : p3 = gmul(p3, p4);
1648 : p4 = gsub(p2, p1);
1649 : } else {
1650 : /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1651 : * I = xr*yi + xi*yr */
1652 : p1 = gmul(p3,p4);
1653 : p3 = gmul(xr,yi);
1654 : p4 = gmul(xi,yr);
1655 : p2 = gsub(p3, p4);
1656 : }
1657 : } else {
1658 : if (did_add(yr, yi, &p4)) {
1659 : /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1660 : * I = xr*yi +xi*yr */
1661 : p1 = gmul(p3,p4);
1662 : p3 = gmul(xr,yi);
1663 : p4 = gmul(xi,yr);
1664 : p2 = gsub(p4, p3);
1665 : } else {
1666 : /* R = xr*yr - xi*yi
1667 : * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1668 : p3 = gneg( gmul(p3, p4) );
1669 : p1 = gmul(xr,yr);
1670 : p2 = gmul(xi,yi);
1671 : p4 = gadd(p1, p2);
1672 :
1673 : p2 = gneg(p2);
1674 : }
1675 : }
1676 : tetpil = avma;
1677 : gel(z,1) = gadd(p1,p2);
1678 : gel(z,2) = gadd(p3,p4);
1679 : #else
1680 260243095 : if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1681 : { /* 3M formula */
1682 559356 : p3 = addii(xr,xi);
1683 559356 : p4 = addii(yr,yi);
1684 559356 : p1 = mulii(xr,yr);
1685 559356 : p2 = mulii(xi,yi);
1686 559356 : p3 = mulii(p3,p4);
1687 559356 : p4 = addii(p2,p1);
1688 559356 : tetpil = avma;
1689 559356 : gel(z,1) = subii(p1,p2);
1690 559356 : gel(z,2) = subii(p3,p4);
1691 559356 : if (!signe(gel(z,2)))
1692 113260 : return gc_INT((pari_sp)(z+3), gel(z,1));
1693 : }
1694 : else
1695 : { /* naive 4M formula: avoid all loss of accuracy */
1696 259683739 : p1 = gmul(xr,yr);
1697 259655359 : p2 = gmul(xi,yi);
1698 259647439 : p3 = gmul(xr,yi);
1699 259658364 : p4 = gmul(xi,yr);
1700 259661227 : tetpil = avma;
1701 259661227 : gel(z,1) = gsub(p1,p2);
1702 259542347 : gel(z,2) = gadd(p3,p4);
1703 259547238 : if (isintzero(gel(z,2)))
1704 : {
1705 50575 : cgiv(gel(z,2));
1706 50575 : return gc_upto((pari_sp)(z+3), gel(z,1));
1707 : }
1708 : }
1709 : #endif
1710 :
1711 259934214 : gc_slice_unsafe(av,tetpil, z+1,2); return z;
1712 : }
1713 : /* x,y PADIC */
1714 : static GEN
1715 17505522 : mulpp(GEN x, GEN y) {
1716 17505522 : GEN pd, p = padic_p(x), ux = padic_u(x), uy = padic_u(y);
1717 17505522 : long e = valp(x) + valp(y), dx, dy;
1718 :
1719 17505522 : if (!equalii(p, padic_p(y))) pari_err_OP("*",x,y);
1720 17505448 : if (!signe(ux) || !signe(uy)) return zeropadic(p, e);
1721 16931530 : dx = precp(x); dy = precp(y);
1722 16931530 : if (dx > dy) { pd = padic_pd(y); dx = dy; } else pd = padic_pd(x);
1723 16931530 : retmkpadic(Fp_mul(ux, uy, pd), icopy(p), icopy(pd), e, dx);
1724 : }
1725 : /* x,y QUAD, w^2 = - b w - c, where b = 0 or -1
1726 : * (ux + vx w)(uy + vy w) = ux uy - c vx vy + (ux vy + uy vx - b vx vy) w */
1727 : static GEN
1728 1127 : mulqq(GEN x, GEN y)
1729 : {
1730 1127 : GEN T = gel(x,1), b = gel(T,3), c = gel(T,2);
1731 1127 : GEN ux = gel(x,2), vx = gel(x,3), uy = gel(y,2), vy = gel(y,3);
1732 : GEN z, s, U, V, E, F;
1733 : pari_sp av, av2;
1734 1127 : z = cgetg(4, t_QUAD), gel(z,1) = ZX_copy(T); av = avma;
1735 1127 : U = gmul(ux, uy);
1736 1127 : V = gmul(vx, vy); s = gmul(c, V);
1737 1127 : E = gmul(gadd(ux, vx), gadd(uy, vy));
1738 1127 : F = signe(b)? U: gadd(U, V);
1739 : /* E - F = ux vy + uy vx - b vx vy */
1740 1127 : av2 = avma;
1741 1127 : gel(z,2) = gsub(U, s);
1742 1127 : gel(z,3) = gsub(E, F);
1743 1127 : gc_slice_unsafe(av, av2, z+2, 2); return z;
1744 : }
1745 :
1746 : GEN
1747 15629259 : mulcxI(GEN x)
1748 : {
1749 : GEN z;
1750 15629259 : switch(typ(x))
1751 : {
1752 1006884 : case t_INT:
1753 1006884 : if (!signe(x)) return gen_0;
1754 : case t_REAL: case t_FRAC:
1755 2610973 : return mkcomplex(gen_0, x);
1756 13018055 : case t_COMPLEX:
1757 13018055 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1758 13004634 : z = cgetg(3,t_COMPLEX);
1759 13004857 : gel(z,1) = gneg(gel(x,2));
1760 13005260 : gel(z,2) = gel(x,1); return z;
1761 56 : default:
1762 56 : return gmul(gen_I(), x);
1763 : }
1764 : }
1765 : GEN
1766 2985897 : mulcxmI(GEN x)
1767 : {
1768 : GEN z;
1769 2985897 : switch(typ(x))
1770 : {
1771 546 : case t_INT:
1772 546 : if (!signe(x)) return gen_0;
1773 : case t_REAL: case t_FRAC:
1774 97845 : return mkcomplex(gen_0, gneg(x));
1775 2685078 : case t_COMPLEX:
1776 2685078 : if (isintzero(gel(x,1))) return gel(x,2);
1777 2676672 : z = cgetg(3,t_COMPLEX);
1778 2676651 : gel(z,1) = gel(x,2);
1779 2676651 : gel(z,2) = gneg(gel(x,1)); return z;
1780 202974 : default:
1781 202974 : return gmul(mkcomplex(gen_0, gen_m1), x);
1782 : }
1783 : }
1784 : /* x * I^k */
1785 : GEN
1786 5563071 : mulcxpowIs(GEN x, long k)
1787 : {
1788 5563071 : switch (k & 3)
1789 : {
1790 1369158 : case 1: return mulcxI(x);
1791 1362214 : case 2: return gneg(x);
1792 1403384 : case 3: return mulcxmI(x);
1793 : }
1794 1428315 : return x;
1795 : }
1796 :
1797 : /* Convert t_POL y to t_SER of length l and valuation e then normalize. Must
1798 : * not be called when y is an exact zero carrying type information, see comment
1799 : * in gmul(t_SER, t_SER): the valuation may incorrectly be altered by
1800 : * normalizeser in this case. Variant of RgX_to_ser_i with a different
1801 : * meaning for e. */
1802 : static GEN
1803 5638655 : pol2ser(GEN y, long l, long e)
1804 : {
1805 5638655 : long i, ly = lg(y);
1806 5638655 : GEN z = cgetg(l, t_SER);
1807 5638655 : z[1] = evalvalser(e) | evalvarn(varn(y)) | evalsigne(1);
1808 5638655 : if (ly >= l) {
1809 25511387 : for (i = 2; i < l; i++) gel(z,i) = gel(y,i);
1810 : } else {
1811 339601 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1812 240504 : for ( ; i < l; i++) gel(z,i) = gen_0;
1813 : }
1814 5638655 : return normalizeser(z);
1815 : }
1816 : /* assume typ(x) = t_VEC */
1817 : static int
1818 98 : is_ext_qfr(GEN x)
1819 91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
1820 189 : && typ(gel(x,2)) == t_REAL; }
1821 :
1822 : GEN
1823 8695138702 : gmul(GEN x, GEN y)
1824 : {
1825 : long tx, ty, lx, ly, vx, vy, i, l;
1826 : pari_sp av, tetpil;
1827 : GEN z, p1;
1828 :
1829 8695138702 : if (x == y) return gsqr(x);
1830 7755936966 : tx = typ(x); ty = typ(y);
1831 7755936966 : if (tx == ty) switch(tx)
1832 : {
1833 3552055650 : case t_INT: return mulii(x,y);
1834 1892578435 : case t_REAL: return mulrr(x,y);
1835 1439078 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1836 1439078 : z = cgetg(3,t_INTMOD);
1837 1439077 : if (X==Y || equalii(X,Y))
1838 1439035 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1839 42 : gel(z,1) = gcdii(X,Y);
1840 42 : warn_coercion(X,Y,gel(z,1));
1841 42 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1842 42 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
1843 : }
1844 22837368 : case t_FRAC:
1845 : {
1846 22837368 : GEN x1 = gel(x,1), x2 = gel(x,2);
1847 22837368 : GEN y1 = gel(y,1), y2 = gel(y,2);
1848 22837368 : z=cgetg(3,t_FRAC);
1849 22837368 : p1 = gcdii(x1, y2);
1850 22837368 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1851 22837367 : p1 = gcdii(x2, y1);
1852 22837366 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1853 22837366 : tetpil = avma;
1854 22837366 : gel(z,2) = mulii(x2,y2);
1855 22837365 : gel(z,1) = mulii(x1,y1);
1856 22837366 : fix_frac_if_int_GC(z,tetpil); return z;
1857 : }
1858 269793182 : case t_COMPLEX: return mulcc(x, y);
1859 12511082 : case t_PADIC: return mulpp(x, y);
1860 882 : case t_QUAD:
1861 882 : if (!ZX_equal(gel(x,1), gel(y,1))) pari_err_OP("*",x,y);
1862 882 : return mulqq(x, y);
1863 13690497 : case t_FFELT: return FF_mul(x, y);
1864 10951746 : case t_POLMOD:
1865 10951746 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1866 8276801 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1867 2674945 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1868 131538731 : case t_POL:
1869 131538731 : vx = varn(x);
1870 131538731 : vy = varn(y);
1871 131538731 : if (vx != vy) {
1872 46850416 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1873 27724518 : else return RgX_Rg_mul(y, x);
1874 : }
1875 84688315 : return RgX_mul(x, y);
1876 :
1877 4431344 : case t_SER:
1878 : {
1879 : long e;
1880 4431344 : vx = varn(x);
1881 4431344 : vy = varn(y);
1882 4431344 : if (vx != vy) {
1883 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1884 3675 : return mul_ser_scal(y, x);
1885 : }
1886 4427669 : e = valser(x) + valser(y);
1887 4427669 : lx = minss(lg(x), lg(y));
1888 4427669 : if (lx == 2) return zeroser(vx, e);
1889 : /* Must handle the case of exact zeros here. If x or y is t^e (z+O(t))
1890 : * where z is an exact zero giving type information, then its valser is
1891 : * v+1 and their product has the same shape but we mustn't call
1892 : * normalizeser() since the result is already normalized (see "dangerous
1893 : * case" there). It would add 1 to the valuation again.
1894 : *
1895 : * But x,y may be zero divisors, Mod(2,4) + O(t) say, and the
1896 : * product again has the above shape but this time we must normalize.
1897 : * There is no way to distinguish between the two situations in
1898 : * pol2ser() from the product of the polynomial parts only */
1899 4427655 : if (lx == 3)
1900 : {
1901 : long s;
1902 : GEN c;
1903 25473 : z = cgetg(3, t_SER); c = gmul(gel(x,2), gel(y,2));
1904 25473 : s = !gequal0(c);
1905 25473 : if (!s && signe(x) && signe(y) && isexactzero(c)) e++;
1906 25473 : z[1] = evalsigne(s) | evalvarn(vx) | evalvalser(e);
1907 25473 : gel(z,2) = c; return z;
1908 : }
1909 4402182 : av = avma;
1910 4402182 : x = ser2pol_i(x, lx);
1911 4402182 : y = ser2pol_i(y, lx);
1912 4402182 : y = RgXn_mul(x, y, lx-2);
1913 4402182 : return gc_GEN(av, pol2ser(y, lx, e));
1914 : }
1915 42 : case t_VEC:
1916 42 : if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
1917 : /* fall through, handle extended t_QFB */
1918 198 : case t_QFB: return qfbcomp(x,y);
1919 6723062 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1920 4127727 : case t_MAT: return RgM_mul(x, y);
1921 :
1922 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1923 1631 : z = cgetg_copy(x, &l);
1924 1631 : if (l != lg(y)) break;
1925 17325 : for (i=1; i<l; i++)
1926 : {
1927 15694 : long yi = y[i];
1928 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1929 15694 : z[i] = x[yi];
1930 : }
1931 1631 : return z;
1932 :
1933 0 : default:
1934 0 : pari_err_TYPE2("*",x,y);
1935 : }
1936 : /* tx != ty */
1937 1835222318 : if (is_const_t(ty) && is_const_t(tx)) {
1938 1672394708 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1939 1672394708 : switch(tx) {
1940 1543110984 : case t_INT:
1941 : switch(ty)
1942 : {
1943 999153225 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1944 1349376 : case t_INTMOD:
1945 1349376 : z = cgetg(3, t_INTMOD);
1946 1349377 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1947 102219744 : case t_FRAC:
1948 102219744 : if (!signe(x)) return gen_0;
1949 85134442 : z=cgetg(3,t_FRAC);
1950 85134611 : p1 = gcdii(x,gel(y,2));
1951 85133850 : if (equali1(p1))
1952 : {
1953 46043405 : set_avma((pari_sp)z);
1954 46043389 : gel(z,2) = icopy(gel(y,2));
1955 46043460 : gel(z,1) = mulii(gel(y,1), x);
1956 : }
1957 : else
1958 : {
1959 39090723 : x = diviiexact(x,p1); tetpil = avma;
1960 39090223 : gel(z,2) = diviiexact(gel(y,2), p1);
1961 39090052 : gel(z,1) = mulii(gel(y,1), x);
1962 39090551 : fix_frac_if_int_GC(z,tetpil);
1963 : }
1964 85134347 : return z;
1965 437356847 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1966 4199397 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1967 1904 : case t_QUAD: return mulRq(x,y);
1968 1607060 : case t_FFELT: return FF_Z_mul(y,x);
1969 : }
1970 :
1971 : case t_REAL:
1972 121477630 : switch(ty)
1973 : {
1974 37740790 : case t_FRAC: return mulrfrac(x, y);
1975 83736798 : case t_COMPLEX: return mulRc(x, y);
1976 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1977 21 : default: pari_err_TYPE2("*",x,y);
1978 : }
1979 :
1980 8288 : case t_INTMOD:
1981 : switch(ty)
1982 : {
1983 7133 : case t_FRAC: { GEN X = gel(x,1);
1984 7133 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1985 7133 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1986 : }
1987 49 : case t_COMPLEX: return mulRc(x,y);
1988 511 : case t_PADIC: { GEN X = gel(x,1);
1989 511 : z = cgetg(3, t_INTMOD);
1990 511 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1991 : }
1992 63 : case t_QUAD: return mulRq(x, y);
1993 532 : case t_FFELT:
1994 532 : if (!equalii(gel(x,1),FF_p_i(y)))
1995 0 : pari_err_OP("*",x,y);
1996 532 : return FF_Z_mul(y,gel(x,2));
1997 : }
1998 :
1999 : case t_FRAC:
2000 : switch(ty)
2001 : {
2002 5372488 : case t_COMPLEX: return mulRc(x, y);
2003 2582761 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
2004 343 : case t_QUAD: return mulRq(x, y);
2005 2079 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
2006 : }
2007 :
2008 : case t_FFELT:
2009 35 : pari_err_TYPE2("*",x,y);
2010 :
2011 21 : case t_COMPLEX:
2012 : switch(ty)
2013 : {
2014 14 : case t_PADIC:
2015 14 : return Zp_nosquare_m1(padic_p(y))? mulRc(y, x): mulTp(x, y);
2016 7 : case t_QUAD:
2017 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
2018 7 : return mulqf(y, x, lx);
2019 : }
2020 :
2021 : case t_PADIC: /* ty == t_QUAD */
2022 28 : return (kro_quad(y,padic_p(x))== -1)? mulRq(x, y): mulTp(y, x);
2023 : }
2024 : }
2025 :
2026 165274898 : if (is_matvec_t(ty)) switch(tx)
2027 : {
2028 121957 : case t_VEC:
2029 : switch(ty) {
2030 14854 : case t_COL: return RgV_RgC_mul(x,y);
2031 107103 : case t_MAT: return RgV_RgM_mul(x,y);
2032 : }
2033 0 : break;
2034 1687 : case t_COL:
2035 : switch(ty) {
2036 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2037 0 : case t_MAT: return RgC_RgM_mul(x,y);
2038 : }
2039 0 : break;
2040 882513 : case t_MAT:
2041 : switch(ty) {
2042 0 : case t_VEC: return RgM_RgV_mul(x,y);
2043 882512 : case t_COL: return RgM_RgC_mul(x,y);
2044 : }
2045 : default:
2046 5517917 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
2047 5517920 : return mul_self_scal(y, x);
2048 : }
2049 161853735 : if (is_matvec_t(tx))
2050 : {
2051 2127852 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2052 2127995 : return mul_self_scal(x, y);
2053 : }
2054 159726029 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2055 : /* tx < ty, !ismatvec(x and y) */
2056 :
2057 159726029 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2058 3215898 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2059 156510131 : if (is_scalar_t(tx)) {
2060 153042091 : if (tx == t_POLMOD) {
2061 3159530 : vx = varn(gel(x,1));
2062 3159530 : vy = gvar(y);
2063 3159530 : if (vx != vy) {
2064 2907313 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2065 23093 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2066 : }
2067 : /* error if ty == t_SER */
2068 252217 : av = avma; y = gmod(y, gel(x,1));
2069 252210 : return gc_upto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2070 : }
2071 149882561 : return mul_scal(y, x, ty);
2072 : }
2073 :
2074 : /* x and y are not scalars, nor matvec */
2075 3468002 : vx = gvar(x);
2076 3468028 : vy = gvar(y);
2077 3468028 : if (vx != vy) /* x or y is treated as a scalar */
2078 2785619 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2079 2785619 : : mul_scal(y, x, ty);
2080 : /* vx = vy */
2081 2044604 : switch(tx)
2082 : {
2083 2044569 : case t_POL:
2084 : switch (ty)
2085 : {
2086 6762 : case t_SER:
2087 : {
2088 : long e, E;
2089 6762 : av = avma; e = RgX_valrem(x, &x);
2090 6762 : if (e == LONG_MAX)
2091 21 : return gc_upto(av, gmul(Rg_get_0(x), Rg_get_0(y)));
2092 6741 : ly = lg(y); E = e + valser(y);
2093 6741 : if (ly == 2) { set_avma(av); return zeroser(vx, E); }
2094 6510 : if (ly == 3)
2095 : {
2096 : long s;
2097 : GEN c;
2098 63 : z = cgetg(3, t_SER); c = gmul(gel(x,2), gel(y,2));
2099 63 : s = !gequal0(c);
2100 63 : if (!s && signe(x) && signe(y) && isexactzero(c)) E++;
2101 63 : z[1] = evalsigne(s) | evalvarn(vx) | evalvalser(E);
2102 63 : gel(z,2) = c; return gc_upto(av, z);
2103 : }
2104 6447 : if (degpol(x))
2105 : {
2106 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2107 2030 : return gc_GEN(av, pol2ser(x, ly, E));
2108 : }
2109 : /* take advantage of x = c*t^e */
2110 4417 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2111 4417 : setvalser(y, e + valser(y)); return y;
2112 : }
2113 :
2114 2037807 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2115 : }
2116 0 : break;
2117 :
2118 35 : case t_SER:
2119 : switch (ty)
2120 : {
2121 35 : case t_RFRAC:
2122 35 : av = avma;
2123 35 : return gc_upto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2124 : }
2125 0 : break;
2126 : }
2127 0 : pari_err_TYPE2("*",x,y);
2128 : return NULL; /* LCOV_EXCL_LINE */
2129 : }
2130 :
2131 : /* return a nonnormalized result */
2132 : GEN
2133 113207 : sqr_ser_part(GEN x, long l1, long l2)
2134 : {
2135 : long i, j, l;
2136 : pari_sp av;
2137 : GEN Z, z, p1, p2;
2138 : long mi;
2139 113207 : if (l2 < l1) return zeroser(varn(x), 2*valser(x));
2140 113193 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2141 113193 : Z = cgetg(l2-l1+3, t_SER);
2142 113193 : Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
2143 113193 : z = Z + 2-l1;
2144 113193 : x += 2; mi = 0;
2145 425866 : for (i=0; i<l1; i++)
2146 : {
2147 312673 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2148 : }
2149 :
2150 721710 : for (i=l1; i<=l2; i++)
2151 : {
2152 608517 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2153 608517 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2154 1257891 : for (j=i-mi; j<=minss(l,mi); j++)
2155 649374 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2156 608517 : p1 = gshift(p1,1);
2157 608517 : if ((i&1) == 0 && p2[i>>1])
2158 94208 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2159 608517 : gel(z,i) = gc_upto(av,p1);
2160 : }
2161 113193 : return Z;
2162 : }
2163 :
2164 : /* (u + v X)^2 mod (X^2 + bX + c), b = 0 or -1
2165 : * = u^2 - c v^2 + (2uv - b v^2) X */
2166 : static GEN
2167 70 : sqrq(GEN x)
2168 : {
2169 70 : GEN T = gel(x,1), c = gel(T,2), b = gel(T,3);
2170 70 : GEN u = gel(x,2), v = gel(x,3), U, V, E, F, s, z;
2171 : pari_sp av, av2;
2172 :
2173 70 : z = cgetg(4, t_QUAD), gel(z,1) = ZX_copy(T); av = avma;
2174 70 : U = gsqr(u);
2175 70 : V = gsqr(v); s = gmul(c, V);
2176 70 : E = gmul(u, v);
2177 70 : F = signe(b)? gadd(E, V): E; /* u v - b v^2 */
2178 70 : av2 = avma;
2179 70 : gel(z,2) = gsub(U, s);
2180 70 : gel(z,3) = gadd(E, F);
2181 70 : gc_slice_unsafe(av, av2, z+2, 2); return z;
2182 : }
2183 :
2184 : GEN
2185 1184033437 : gsqr(GEN x)
2186 : {
2187 : long i, lx;
2188 : pari_sp av, tetpil;
2189 : GEN z, p1, p2, p3;
2190 :
2191 1184033437 : switch(typ(x))
2192 : {
2193 971236093 : case t_INT: return sqri(x);
2194 196314015 : case t_REAL: return sqrr(x);
2195 142500 : case t_INTMOD: { GEN X = gel(x,1);
2196 142500 : z = cgetg(3,t_INTMOD);
2197 142500 : gel(z,2) = gc_INT((pari_sp)z, remii(sqri(gel(x,2)), X));
2198 142497 : gel(z,1) = icopy(X); return z;
2199 : }
2200 4146871 : case t_FRAC: return sqrfrac(x);
2201 :
2202 7505913 : case t_COMPLEX:
2203 7505913 : if (isintzero(gel(x,1))) {
2204 239995 : av = avma;
2205 239995 : return gc_upto(av, gneg(gsqr(gel(x,2))));
2206 : }
2207 7265915 : z = cgetg(3,t_COMPLEX); av = avma;
2208 7265898 : p1 = gadd(gel(x,1),gel(x,2));
2209 7265742 : p2 = gsub(gel(x,1), gel(x,2));
2210 7265762 : p3 = gmul(gel(x,1),gel(x,2));
2211 7265875 : tetpil = avma;
2212 7265875 : gel(z,1) = gmul(p1,p2);
2213 7265920 : gel(z,2) = gshift(p3,1); gc_slice_unsafe(av,tetpil,z+1,2); return z;
2214 :
2215 4795 : case t_PADIC:
2216 : {
2217 4795 : GEN u = padic_u(x), p = padic_p(x), pd = padic_pd(x);
2218 4795 : long v2 = 2*valp(x), d = precp(x);
2219 4795 : if (!signe(u)) { x = gcopy(x); setvalp(x, v2); return x; }
2220 4739 : if (!absequaliu(p, 2))
2221 3206 : retmkpadic(Fp_sqr(u, pd), icopy(p), icopy(pd), v2, d);
2222 : /* p = 2*/
2223 1533 : if (d == 1) /* (1 + O(2))^2 = 1 + O(2^3) */
2224 7 : retmkpadic(gen_1, gen_2, utoipos(8), v2, 3);
2225 1526 : retmkpadic(remi2n(sqri(u), d + 1), gen_2, int2n(d + 1), v2, d + 1);
2226 : }
2227 70 : case t_QUAD: return sqrq(x);
2228 :
2229 419385 : case t_POLMOD: return sqr_polmod(gel(x,1), gel(x,2));
2230 :
2231 2970752 : case t_FFELT: return FF_sqr(x);
2232 :
2233 1279174 : case t_POL: return RgX_sqr(x);
2234 :
2235 35490 : case t_SER:
2236 35490 : lx = lg(x);
2237 35490 : if (ser_isexactzero(x)) {
2238 21 : GEN z = gcopy(x);
2239 21 : setvalser(z, 2*valser(x));
2240 21 : return z;
2241 : }
2242 35469 : if (lx < 40)
2243 35210 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2244 : else
2245 : {
2246 259 : pari_sp av = avma;
2247 259 : long v = 2 * valser(x);
2248 259 : x = ser2pol_i(x, lx);
2249 259 : x = RgXn_sqr(x, lx-2);
2250 259 : return gc_GEN(av, pol2ser(x, lx, v));
2251 : }
2252 :
2253 63 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2254 63 : gel(z,1) = gsqr(gel(x,1));
2255 63 : gel(z,2) = gsqr(gel(x,2)); return z;
2256 :
2257 1169 : case t_MAT: return RgM_sqr(x);
2258 28 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
2259 : /* fall through handle extended t_QFB */
2260 36 : case t_QFB: return qfbsqr(x);
2261 658 : case t_VECSMALL:
2262 658 : z = cgetg_copy(x, &lx);
2263 16289 : for (i=1; i<lx; i++)
2264 : {
2265 15631 : long xi = x[i];
2266 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2267 15631 : z[i] = x[xi];
2268 : }
2269 658 : return z;
2270 : }
2271 7 : pari_err_TYPE2("*",x,x);
2272 : return NULL; /* LCOV_EXCL_LINE */
2273 : }
2274 :
2275 : /********************************************************************/
2276 : /** **/
2277 : /** DIVISION **/
2278 : /** **/
2279 : /********************************************************************/
2280 : static GEN
2281 32139 : div_rfrac_scal(GEN x, GEN y)
2282 : {
2283 32139 : pari_sp av = avma;
2284 32139 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2285 32139 : return gc_upto(av, gred_rfrac_simple(gel(x,1), d));
2286 : }
2287 : static GEN
2288 37767 : div_scal_rfrac(GEN x, GEN y)
2289 : {
2290 37767 : GEN y1 = gel(y,1), y2 = gel(y,2);
2291 37767 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2292 : {
2293 14 : if (degpol(y1))
2294 : {
2295 14 : pari_sp av = avma;
2296 14 : GEN _1 = Rg_get_1(x);
2297 14 : if (transtype(_1)) y1 = gmul(y1, _1);
2298 14 : return gc_upto(av, gred_rfrac_simple(gmul(x, y2), y1));
2299 : }
2300 0 : y1 = gel(y1,2);
2301 : }
2302 37753 : return RgX_Rg_mul(y2, gdiv(x,y1));
2303 : }
2304 : static GEN
2305 1187237 : div_rfrac(GEN x, GEN y)
2306 1187237 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2307 :
2308 : /* x != 0 */
2309 : static GEN
2310 1343391 : div_ser_scal(GEN x, GEN t)
2311 : {
2312 1343391 : if (ser_isexactzero(x))
2313 : {
2314 28 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gdiv(gel(x,2), t), varn(x), 1);
2315 28 : setvalser(z, valser(x)); return z;
2316 : }
2317 6223541 : pari_APPLY_ser(gdiv(gel(x,i), t));
2318 : }
2319 : GEN
2320 7658 : ser_normalize(GEN x)
2321 : {
2322 7658 : long i, lx = lg(x);
2323 : GEN c, z;
2324 7658 : if (lx == 2) return x;
2325 7658 : c = gel(x,2); if (gequal1(c)) return x;
2326 7581 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2327 108752 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2328 7581 : return z;
2329 : }
2330 :
2331 : /* y != 0 */
2332 : static GEN
2333 7128718 : div_T_scal(GEN x, GEN y, long tx) {
2334 7128718 : switch(tx)
2335 : {
2336 5757039 : case t_POL: return RgX_Rg_div(x, y);
2337 1343384 : case t_SER: return div_ser_scal(x, y);
2338 28296 : case t_RFRAC: return div_rfrac_scal(x,y);
2339 : }
2340 0 : pari_err_TYPE2("/",x,y);
2341 : return NULL; /* LCOV_EXCL_LINE */
2342 : }
2343 :
2344 : static GEN
2345 9258206 : div_scal_pol(GEN x, GEN y) {
2346 9258206 : long ly = lg(y);
2347 : pari_sp av;
2348 : GEN _1;
2349 9258206 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2350 9179861 : if (isrationalzero(x)) return zeropol(varn(y));
2351 7133077 : av = avma;
2352 7133077 : _1 = Rg_get_1(x); if (transtype(_1)) y = gmul(y, _1);
2353 7133077 : return gc_upto(av, gred_rfrac_simple(x,y));
2354 : }
2355 : static GEN
2356 19096 : div_scal_ser(GEN x, GEN y)
2357 : {
2358 19096 : pari_sp av = avma;
2359 19096 : GEN _1 = Rg_get_1(x);
2360 19096 : if (transtype(_1)) y = gmul(y, _1);
2361 19096 : return gc_upto(av, gmul(x, ser_inv(y)));
2362 : }
2363 : static GEN
2364 9266449 : div_scal_T(GEN x, GEN y, long ty) {
2365 9266449 : switch(ty)
2366 : {
2367 9210349 : case t_POL: return div_scal_pol(x, y);
2368 19089 : case t_SER: return div_scal_ser(x, y);
2369 37011 : case t_RFRAC: return div_scal_rfrac(x, y);
2370 : }
2371 0 : pari_err_TYPE2("/",x,y);
2372 : return NULL; /* LCOV_EXCL_LINE */
2373 : }
2374 :
2375 : /* assume tx = ty = t_SER, same variable vx */
2376 : static GEN
2377 774266 : div_ser(GEN x, GEN y, long vx)
2378 : {
2379 774266 : long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
2380 774266 : GEN y0 = y, z;
2381 774266 : pari_sp av = avma;
2382 :
2383 774266 : if (!signe(y)) pari_err_INV("div_ser", y);
2384 774259 : if (ser_isexactzero(x))
2385 : {
2386 59892 : if (lx == 2) return zeroser(vx, v);
2387 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2388 147 : setvalser(z, v); return z;
2389 : }
2390 714367 : if (lx < ly) ly = lx;
2391 714367 : y = ser2pol_i_normalize(y, ly, &e);
2392 714367 : if (e)
2393 : {
2394 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2395 7 : ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
2396 : }
2397 714367 : if (ly == 3)
2398 : {
2399 12452 : z = cgetg(3, t_SER);
2400 12452 : z[1] = evalsigne(1) | evalvalser(v) | evalvarn(vx);
2401 12452 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2402 12452 : if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
2403 12452 : return gc_upto(av, z);
2404 : }
2405 701915 : x = ser2pol_i(x, ly);
2406 701915 : y = RgXn_div_i(x, y, ly-2);
2407 701915 : return gc_GEN(av, pol2ser(y, ly, v));
2408 : }
2409 : /* x,y compatible PADIC */
2410 : static GEN
2411 13212814 : divpp(GEN x, GEN y)
2412 : {
2413 13212814 : GEN M, ux = padic_u(x), uy = padic_u(y), p = padic_p(x);
2414 13212814 : long a = precp(x), b = precp(y), v = valp(x) - valp(y);
2415 :
2416 13212814 : if (!signe(uy)) pari_err_INV("divpp",y);
2417 13212835 : if (!signe(ux)) return zeropadic(p, v);
2418 13212506 : if (a > b) M = padic_pd(y); else { M = padic_pd(x); b = a; }
2419 13212506 : retmkpadic(Fp_div(ux, uy, M), icopy(padic_p(x)), icopy(M), v, b);
2420 : }
2421 : static GEN
2422 37016 : div_polmod_same(GEN T, GEN x, GEN y)
2423 : {
2424 37016 : long v = varn(T);
2425 37016 : GEN a, z = cgetg(3, t_POLMOD);
2426 37016 : gel(z,1) = RgX_copy(T);
2427 37016 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2428 11354 : a = gdiv(x, y);
2429 25662 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2430 105 : {
2431 105 : pari_sp av = avma;
2432 105 : a = gc_upto(av, gmul(x, RgXQ_inv(y, T)));
2433 : }
2434 25557 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2435 5278 : {
2436 5278 : pari_sp av = avma;
2437 5278 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2438 5278 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2439 5278 : a = gc_upto(av, a);
2440 : }
2441 : else
2442 : {
2443 20279 : pari_sp av = avma;
2444 20279 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2445 20279 : a = gc_upto(av, a);
2446 : }
2447 37016 : gel(z,2) = a; return z;
2448 : }
2449 : GEN
2450 439096338 : gdiv(GEN x, GEN y)
2451 : {
2452 439096338 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2453 : pari_sp av, tetpil;
2454 : GEN z, p1;
2455 :
2456 439096338 : if (tx == ty) switch(tx)
2457 : {
2458 92088877 : case t_INT:
2459 92088877 : return Qdivii(x,y);
2460 :
2461 94285714 : case t_REAL: return divrr(x,y);
2462 25277 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2463 25277 : z = cgetg(3,t_INTMOD);
2464 25277 : if (X==Y || equalii(X,Y))
2465 25263 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2466 14 : gel(z,1) = gcdii(X,Y);
2467 14 : warn_coercion(X,Y,gel(z,1));
2468 14 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2469 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
2470 : }
2471 4044157 : case t_FRAC: {
2472 4044157 : GEN x1 = gel(x,1), x2 = gel(x,2);
2473 4044157 : GEN y1 = gel(y,1), y2 = gel(y,2);
2474 4044157 : z = cgetg(3, t_FRAC);
2475 4044157 : p1 = gcdii(x1, y1);
2476 4044155 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2477 4044155 : p1 = gcdii(x2, y2);
2478 4044156 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2479 4044151 : tetpil = avma;
2480 4044151 : gel(z,2) = mulii(x2,y1);
2481 4044155 : gel(z,1) = mulii(x1,y2);
2482 4044154 : normalize_frac(z);
2483 4044154 : fix_frac_if_int_GC(z,tetpil);
2484 4044157 : return z;
2485 : }
2486 9426628 : case t_COMPLEX:
2487 9426628 : if (isintzero(gel(y,1)))
2488 : {
2489 182777 : y = gel(y,2);
2490 182777 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2491 103943 : z = cgetg(3,t_COMPLEX);
2492 103943 : gel(z,1) = gdiv(gel(x,2), y);
2493 103943 : av = avma;
2494 103943 : gel(z,2) = gc_upto(av, gneg(gdiv(gel(x,1), y)));
2495 103943 : return z;
2496 : }
2497 9243845 : av = avma;
2498 9243845 : return gc_upto(av, gdiv(mulcc(x, conj_i(y)), cxnorm(y)));
2499 :
2500 587462 : case t_PADIC:
2501 587462 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("/",x,y);
2502 587462 : return divpp(x, y);
2503 :
2504 259 : case t_QUAD:
2505 259 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2506 245 : av = avma;
2507 245 : return gc_upto(av, gdiv(mulqq(x, conj_i(y)), quadnorm(y)));
2508 :
2509 246148 : case t_FFELT: return FF_div(x,y);
2510 :
2511 41391 : case t_POLMOD:
2512 41391 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2513 37016 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2514 : else {
2515 4375 : av = avma;
2516 4375 : z = gc_upto(av, gmul(x, ginv(y)));
2517 : }
2518 41391 : return z;
2519 :
2520 18569301 : case t_POL:
2521 18569301 : vx = varn(x);
2522 18569301 : vy = varn(y);
2523 18569301 : if (vx != vy) {
2524 154175 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2525 47857 : else return div_scal_pol(x, y);
2526 : }
2527 18415126 : if (!signe(y)) pari_err_INV("gdiv",y);
2528 18415126 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2529 18287066 : av = avma;
2530 18287066 : return gc_upto(av, gred_rfrac2(x,y));
2531 :
2532 774287 : case t_SER:
2533 774287 : vx = varn(x);
2534 774287 : vy = varn(y);
2535 774287 : if (vx != vy) {
2536 21 : if (varncmp(vx, vy) < 0)
2537 : {
2538 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2539 7 : return div_ser_scal(x, y);
2540 : }
2541 7 : return div_scal_ser(x, y);
2542 : }
2543 774266 : return div_ser(x, y, vx);
2544 1191696 : case t_RFRAC:
2545 1191696 : vx = varn(gel(x,2));
2546 1191696 : vy = varn(gel(y,2));
2547 1191696 : if (vx != vy) {
2548 4459 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2549 756 : else return div_scal_rfrac(x, y);
2550 : }
2551 1187237 : return div_rfrac(x,y);
2552 :
2553 21 : case t_VEC: /* handle extended t_QFB */
2554 21 : case t_QFB: av = avma; return gc_upto(av, qfbcomp(x, ginv(y)));
2555 :
2556 28 : case t_MAT:
2557 28 : p1 = RgM_div(x,y);
2558 28 : if (!p1) pari_err_INV("gdiv",y);
2559 21 : return p1;
2560 :
2561 0 : default: pari_err_TYPE2("/",x,y);
2562 : }
2563 :
2564 217818533 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2565 : {
2566 3808834 : long s = signe(x);
2567 3808834 : if (!s) {
2568 522133 : if (gequal0(y)) pari_err_INV("gdiv",y);
2569 522133 : switch (ty)
2570 : {
2571 518990 : default: return gen_0;
2572 28 : case t_INTMOD:
2573 28 : z = cgetg(3,t_INTMOD);
2574 28 : gel(z,1) = icopy(gel(y,1));
2575 28 : gel(z,2) = gen_0; return z;
2576 3115 : case t_FFELT: return FF_zero(y);
2577 : }
2578 : }
2579 3286701 : if (is_pm1(x)) {
2580 1323860 : if (s > 0) return ginv(y);
2581 228936 : av = avma; return gc_upto(av, ginv(gneg(y)));
2582 : }
2583 1962841 : switch(ty)
2584 : {
2585 646520 : case t_REAL: return divir(x,y);
2586 21 : case t_INTMOD:
2587 21 : z = cgetg(3, t_INTMOD);
2588 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2589 792522 : case t_FRAC:
2590 792522 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2591 792522 : if (equali1(p1))
2592 : {
2593 271596 : set_avma((pari_sp)z);
2594 271596 : gel(z,2) = icopy(gel(y,1));
2595 271596 : gel(z,1) = mulii(gel(y,2), x);
2596 271596 : normalize_frac(z);
2597 271596 : fix_frac_if_int(z);
2598 : }
2599 : else
2600 : {
2601 520926 : x = diviiexact(x,p1); tetpil = avma;
2602 520926 : gel(z,2) = diviiexact(gel(y,1), p1);
2603 520926 : gel(z,1) = mulii(gel(y,2), x);
2604 520926 : normalize_frac(z);
2605 520926 : fix_frac_if_int_GC(z,tetpil);
2606 : }
2607 792522 : return z;
2608 :
2609 420 : case t_FFELT: return Z_FF_div(x,y);
2610 521629 : case t_COMPLEX: return divRc(x,y);
2611 1505 : case t_PADIC: return divTp(x, y);
2612 231 : case t_QUAD: return divRq(x,y);
2613 : }
2614 : }
2615 214009691 : if (gequal0(y))
2616 : {
2617 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2618 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2619 : }
2620 :
2621 214016510 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2622 : {
2623 173391520 : case t_REAL:
2624 173391520 : switch(ty)
2625 : {
2626 171108008 : case t_INT: return divri(x,y);
2627 519570 : case t_FRAC:
2628 519570 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2629 519570 : return gc_leaf(av, z);
2630 1763911 : case t_COMPLEX: return divRc(x, y);
2631 42 : case t_QUAD: return divfq(x, y, realprec(x));
2632 18 : default: pari_err_TYPE2("/",x,y);
2633 : }
2634 :
2635 294 : case t_INTMOD:
2636 : switch(ty)
2637 : {
2638 210 : case t_INT:
2639 210 : z = cgetg(3, t_INTMOD);
2640 210 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2641 28 : case t_FRAC: { GEN X = gel(x,1);
2642 28 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2643 28 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2644 : }
2645 7 : case t_FFELT:
2646 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2647 0 : pari_err_OP("/",x,y);
2648 7 : return Z_FF_div(gel(x,2),y);
2649 :
2650 0 : case t_COMPLEX: return divRc(x,y);
2651 14 : case t_QUAD: return divRq(x,y);
2652 :
2653 7 : case t_PADIC: { GEN X = gel(x,1);
2654 7 : z = cgetg(3, t_INTMOD);
2655 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2656 : }
2657 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2658 : }
2659 :
2660 : case t_FRAC:
2661 : switch(ty)
2662 : {
2663 2102207 : case t_INT: z = cgetg(3, t_FRAC);
2664 2102207 : p1 = gcdii(y,gel(x,1));
2665 2102207 : if (equali1(p1))
2666 : {
2667 834898 : set_avma((pari_sp)z); tetpil = 0;
2668 834898 : gel(z,1) = icopy(gel(x,1));
2669 : }
2670 : else
2671 : {
2672 1267309 : y = diviiexact(y,p1); tetpil = avma;
2673 1267309 : gel(z,1) = diviiexact(gel(x,1), p1);
2674 : }
2675 2102207 : gel(z,2) = mulii(gel(x,2),y);
2676 2102207 : normalize_frac(z);
2677 2102207 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2678 2102207 : return z;
2679 :
2680 59572 : case t_REAL:
2681 59572 : av = avma;
2682 59572 : return gc_leaf(av, divir(gel(x,1), mulri(y,gel(x,2))));
2683 :
2684 7 : case t_INTMOD: { GEN Y = gel(y,1);
2685 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2686 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2687 : }
2688 :
2689 28 : case t_FFELT: av=avma;
2690 28 : return gc_upto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2691 :
2692 868 : case t_COMPLEX: return divRc(x, y);
2693 :
2694 2141 : case t_PADIC:
2695 2141 : if (!signe(gel(x,1))) return gen_0;
2696 2141 : return divTp(x, y);
2697 :
2698 42 : case t_QUAD: return divRq(x, y);
2699 : }
2700 :
2701 : case t_FFELT:
2702 161 : switch (ty)
2703 : {
2704 77 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2705 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2706 7 : case t_INTMOD:
2707 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2708 0 : pari_err_OP("/",x,y);
2709 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2710 49 : default:
2711 49 : pari_err_TYPE2("/",x,y);
2712 : }
2713 0 : break;
2714 :
2715 13735528 : case t_COMPLEX:
2716 : switch(ty)
2717 : {
2718 13735532 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2719 0 : case t_INTMOD: return mulRc(ginv(y), x);
2720 0 : case t_PADIC:
2721 0 : return Zp_nosquare_m1(padic_p(y))? divcR(x,y): divTp(x, y);
2722 0 : case t_QUAD:
2723 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2724 0 : return divfq(x, y, lx);
2725 : }
2726 :
2727 : case t_PADIC:
2728 : switch(ty)
2729 : {
2730 1205001 : case t_INT: case t_FRAC: { GEN p = padic_p(x);
2731 1204910 : return signe(padic_u(x))? divpT(x, y)
2732 2409911 : : zeropadic(p, valp(x) - Q_pval(y,p));
2733 : }
2734 7 : case t_INTMOD: { GEN Y = gel(y,1);
2735 7 : z = cgetg(3, t_INTMOD);
2736 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2737 : }
2738 0 : case t_COMPLEX: return divRc(x,y);
2739 14 : case t_QUAD: return divRq(x,y);
2740 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2741 : }
2742 :
2743 : case t_QUAD:
2744 : switch (ty)
2745 : {
2746 1204 : case t_INT: case t_INTMOD: case t_FRAC:
2747 1204 : z = cgetg(4,t_QUAD);
2748 1204 : gel(z,1) = ZX_copy(gel(x,1));
2749 1204 : gel(z,2) = gdiv(gel(x,2), y);
2750 1204 : gel(z,3) = gdiv(gel(x,3), y); return z;
2751 63 : case t_REAL: return divqf(x, y, realprec(y));
2752 28 : case t_PADIC: return divTp(x, y);
2753 0 : case t_COMPLEX:
2754 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2755 0 : return divqf(x, y, ly);
2756 : }
2757 : }
2758 23517575 : switch(ty) {
2759 584932 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2760 584932 : av = avma; return gc_upto(av, gmul(x, ginv(y)));
2761 42 : case t_MAT:
2762 42 : av = avma; p1 = RgM_inv(y);
2763 35 : if (!p1) pari_err_INV("gdiv",y);
2764 28 : return gc_upto(av, gmul(x, p1));
2765 0 : case t_VEC: case t_COL:
2766 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2767 0 : pari_err_TYPE2("/",x,y);
2768 : }
2769 22932601 : switch(tx) {
2770 4427591 : case t_VEC: case t_COL: case t_MAT:
2771 14382579 : pari_APPLY_same(gdiv(gel(x,i),y));
2772 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2773 0 : pari_err_TYPE2("/",x,y);
2774 : }
2775 :
2776 18505010 : vy = gvar(y);
2777 18505769 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2778 218121 : vx = varn(X);
2779 218121 : if (vx != vy) {
2780 217351 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2781 216938 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2782 : }
2783 : /* y is POL, SER or RFRAC */
2784 770 : av = avma;
2785 770 : switch(ty)
2786 : {
2787 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2788 770 : default: y = ginvmod(gmod(y,X), X);
2789 : }
2790 763 : return gc_upto(av, mul_polmod_same(X, gel(x,2), y));
2791 : }
2792 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2793 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2794 : * variable NO_VARIABLE, then the operation is incorrect */
2795 18287648 : vx = gvar(x);
2796 18287641 : if (vx != vy) { /* includes cases where one is scalar */
2797 16394760 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2798 9266036 : else return div_scal_T(x, y, ty);
2799 : }
2800 1892881 : switch(tx)
2801 : {
2802 1046171 : case t_POL:
2803 : switch(ty)
2804 : {
2805 28 : case t_SER:
2806 : {
2807 28 : GEN y0 = y;
2808 : long v;
2809 28 : av = avma; v = RgX_valrem(x, &x);
2810 28 : if (v == LONG_MAX) return gc_upto(av, Rg_get_0(x));
2811 14 : v -= valser(y); ly = lg(y); /* > 2 */
2812 14 : y = ser2pol_i_normalize(y, ly, &i);
2813 14 : if (i)
2814 : {
2815 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2816 7 : ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
2817 : }
2818 7 : x = RgXn_div(x, y, ly-2);
2819 7 : return gc_GEN(av, pol2ser(x, ly, v));
2820 : }
2821 :
2822 1046143 : case t_RFRAC:
2823 : {
2824 1046143 : GEN y1 = gel(y,1), y2 = gel(y,2);
2825 1046143 : if (typ(y1) == t_POL && varn(y1) == vx)
2826 1143 : return mul_rfrac_scal(y2, y1, x);
2827 1045000 : av = avma;
2828 1045000 : return gc_upto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2829 : }
2830 : }
2831 0 : break;
2832 :
2833 584324 : case t_SER:
2834 : switch(ty)
2835 : {
2836 584310 : case t_POL:
2837 : {
2838 584310 : long e = valser(x);
2839 584310 : lx = lg(x);
2840 584310 : if (lx == 2) return zeroser(vx, e - RgX_val(y));
2841 584205 : av = avma; e -= RgX_valrem_inexact(y, &y);
2842 584205 : if (lx == 3)
2843 : {
2844 : GEN c;
2845 51943 : z = cgetg(3, t_SER); c = gdiv(gel(x,2), gel(y,2));
2846 51943 : z[1] = evalsigne(!gequal0(c)) | evalvarn(vx) | evalvalser(e);
2847 51943 : gel(z,2) = c; return z;
2848 : }
2849 532262 : x = ser2pol_i(x, lx);
2850 532262 : x = RgXn_div(x, y, lx - 2);
2851 532262 : return gc_GEN(av, pol2ser(x, lx, e));
2852 : }
2853 14 : case t_RFRAC:
2854 14 : av = avma;
2855 14 : return gc_upto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2856 : }
2857 0 : break;
2858 :
2859 262365 : case t_RFRAC:
2860 : switch(ty)
2861 : {
2862 262365 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2863 0 : case t_SER:
2864 0 : av = avma;
2865 0 : return gc_upto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2866 : }
2867 0 : break;
2868 : }
2869 21 : pari_err_TYPE2("/",x,y);
2870 : return NULL; /* LCOV_EXCL_LINE */
2871 : }
2872 :
2873 : /********************************************************************/
2874 : /** **/
2875 : /** SIMPLE MULTIPLICATION **/
2876 : /** **/
2877 : /********************************************************************/
2878 : GEN
2879 239693647 : gmulsg(long s, GEN x)
2880 : {
2881 : pari_sp av;
2882 : long i;
2883 : GEN z;
2884 :
2885 239693647 : switch(typ(x))
2886 : {
2887 145633558 : case t_INT: return mulsi(s,x);
2888 70308975 : case t_REAL: return s? mulsr(s,x): gen_0; /* gmul semantic */
2889 174874 : case t_INTMOD: { GEN p = gel(x,1);
2890 174874 : z = cgetg(3,t_INTMOD);
2891 174873 : gel(z,2) = gc_INT((pari_sp)z, modii(mulsi(s,gel(x,2)), p));
2892 174870 : gel(z,1) = icopy(p); return z;
2893 : }
2894 548086 : case t_FFELT: return FF_Z_mul(x,stoi(s));
2895 6158995 : case t_FRAC:
2896 6158995 : if (!s) return gen_0;
2897 6148299 : z = cgetg(3,t_FRAC);
2898 6148299 : i = labs(s); i = ugcdiu(gel(x,2), i);
2899 6148299 : if (i == 1)
2900 : {
2901 2259548 : gel(z,2) = icopy(gel(x,2));
2902 2259548 : gel(z,1) = mulis(gel(x,1), s);
2903 : }
2904 : else
2905 : {
2906 3888751 : gel(z,2) = diviuexact(gel(x,2), (ulong)i);
2907 3888751 : gel(z,1) = mulis(gel(x,1), s/i);
2908 3888751 : fix_frac_if_int(z);
2909 : }
2910 6148299 : return z;
2911 :
2912 14542419 : case t_COMPLEX:
2913 14542419 : if (!s) return gen_0;
2914 13416615 : z = cgetg(3, t_COMPLEX);
2915 13416575 : gel(z,1) = gmulsg(s,gel(x,1));
2916 13415079 : gel(z,2) = gmulsg(s,gel(x,2)); return z;
2917 :
2918 1449 : case t_PADIC:
2919 1449 : if (!s) return gen_0;
2920 1449 : av = avma; return gc_upto(av, mulpp(cvtop2(stoi(s),x), x));
2921 :
2922 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2923 7 : gel(z,1) = ZX_copy(gel(x,1));
2924 7 : gel(z,2) = gmulsg(s,gel(x,2));
2925 7 : gel(z,3) = gmulsg(s,gel(x,3)); return z;
2926 :
2927 208516 : case t_POLMOD:
2928 208516 : retmkpolmod(gmulsg(s,gel(x,2)), RgX_copy(gel(x,1)));
2929 :
2930 821141 : case t_POL:
2931 821141 : if (!signe(x)) return RgX_copy(x);
2932 801079 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2933 3327329 : pari_APPLY_pol(gmulsg(s,gel(x,i)));
2934 :
2935 182 : case t_SER:
2936 182 : if (ser_isexactzero(x)) return gcopy(x);
2937 182 : if (!s) return Rg_get_0(x);
2938 3864 : pari_APPLY_ser(gmulsg(s,gel(x,i)));
2939 :
2940 0 : case t_RFRAC:
2941 0 : if (!s) return zeropol(varn(gel(x,2)));
2942 0 : if (s == 1) return gcopy(x);
2943 0 : if (s == -1) return gneg(x);
2944 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), stoi(s));
2945 :
2946 1304779 : case t_VEC: case t_COL: case t_MAT:
2947 4110404 : pari_APPLY_same(gmulsg(s,gel(x,i)));
2948 : }
2949 0 : pari_err_TYPE("gmulsg",x);
2950 : return NULL; /* LCOV_EXCL_LINE */
2951 : }
2952 :
2953 : GEN
2954 283683119 : gmulug(ulong s, GEN x)
2955 : {
2956 : pari_sp av;
2957 : long i;
2958 : GEN z;
2959 :
2960 283683119 : switch(typ(x))
2961 : {
2962 281367155 : case t_INT: return mului(s,x);
2963 1065636 : case t_REAL: return s? mulur(s,x): gen_0; /* gmul semantic */
2964 364 : case t_INTMOD: { GEN p = gel(x,1);
2965 364 : z = cgetg(3,t_INTMOD);
2966 364 : gel(z,2) = gc_INT((pari_sp)z, modii(mului(s,gel(x,2)), p));
2967 364 : gel(z,1) = icopy(p); return z;
2968 : }
2969 413 : case t_FFELT: return FF_Z_mul(x,utoi(s));
2970 981585 : case t_FRAC:
2971 981585 : if (!s) return gen_0;
2972 981571 : z = cgetg(3,t_FRAC);
2973 981571 : i = ugcdiu(gel(x,2), s);
2974 981571 : if (i == 1)
2975 : {
2976 822559 : gel(z,2) = icopy(gel(x,2));
2977 822559 : gel(z,1) = muliu(gel(x,1), s);
2978 : }
2979 : else
2980 : {
2981 159012 : gel(z,2) = diviuexact(gel(x,2), i);
2982 159012 : gel(z,1) = muliu(gel(x,1), s/i);
2983 159012 : fix_frac_if_int(z);
2984 : }
2985 981571 : return z;
2986 :
2987 35945 : case t_COMPLEX:
2988 35945 : if (!s) return gen_0;
2989 35945 : z = cgetg(3, t_COMPLEX);
2990 35945 : gel(z,1) = gmulug(s,gel(x,1));
2991 35945 : gel(z,2) = gmulug(s,gel(x,2)); return z;
2992 :
2993 7342 : case t_PADIC:
2994 7342 : if (!s) return gen_0;
2995 7342 : av = avma; return gc_upto(av, mulpp(cvtop2(utoi(s),x), x));
2996 :
2997 0 : case t_QUAD: z = cgetg(4, t_QUAD);
2998 0 : gel(z,1) = ZX_copy(gel(x,1));
2999 0 : gel(z,2) = gmulug(s,gel(x,2));
3000 0 : gel(z,3) = gmulug(s,gel(x,3)); return z;
3001 :
3002 6363 : case t_POLMOD:
3003 6363 : retmkpolmod(gmulug(s,gel(x,2)), RgX_copy(gel(x,1)));
3004 :
3005 200697 : case t_POL:
3006 200697 : if (!signe(x)) return RgX_copy(x);
3007 199920 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
3008 621229 : pari_APPLY_pol(gmulug(s,gel(x,i)));
3009 :
3010 0 : case t_SER:
3011 0 : if (ser_isexactzero(x)) return gcopy(x);
3012 0 : if (!s) return Rg_get_0(x);
3013 0 : pari_APPLY_ser(gmulug(s,gel(x,i)));
3014 :
3015 0 : case t_RFRAC:
3016 0 : if (!s) return zeropol(varn(gel(x,2)));
3017 0 : if (s == 1) return gcopy(x);
3018 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), utoi(s));
3019 :
3020 17619 : case t_VEC: case t_COL: case t_MAT:
3021 96495 : pari_APPLY_same(gmulug(s,gel(x,i)));
3022 : }
3023 0 : pari_err_TYPE("gmulsg",x);
3024 : return NULL; /* LCOV_EXCL_LINE */
3025 : }
3026 :
3027 : /********************************************************************/
3028 : /** **/
3029 : /** SIMPLE DIVISION **/
3030 : /** **/
3031 : /********************************************************************/
3032 :
3033 : GEN
3034 13001843 : gdivgs(GEN x, long s)
3035 : {
3036 13001843 : long tx = typ(x), i;
3037 : pari_sp av;
3038 : GEN z;
3039 :
3040 13001843 : if (!s)
3041 : {
3042 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3043 0 : pari_err_INV("gdivgs",gen_0);
3044 : }
3045 13001846 : switch(tx)
3046 : {
3047 1592288 : case t_INT: return Qdivis(x, s);
3048 8559614 : case t_REAL: return divrs(x,s);
3049 :
3050 357 : case t_INTMOD:
3051 357 : z = cgetg(3, t_INTMOD);
3052 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3053 :
3054 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3055 :
3056 556880 : case t_FRAC: z = cgetg(3, t_FRAC);
3057 556880 : i = ugcdiu(gel(x,1), labs(s));
3058 556880 : if (i == 1)
3059 : {
3060 413495 : gel(z,2) = mulsi(s, gel(x,2));
3061 413495 : gel(z,1) = icopy(gel(x,1));
3062 : }
3063 : else
3064 : {
3065 143385 : gel(z,2) = mulsi(s/i, gel(x,2));
3066 143385 : gel(z,1) = divis(gel(x,1), i);
3067 : }
3068 556880 : normalize_frac(z);
3069 556880 : fix_frac_if_int(z); return z;
3070 :
3071 1809640 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3072 1809640 : gel(z,1) = gdivgs(gel(x,1),s);
3073 1809640 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3074 :
3075 133 : case t_PADIC: /* divpT */
3076 : {
3077 133 : GEN p = padic_p(x);
3078 133 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3079 133 : av = avma;
3080 133 : return gc_upto(av, divpp(x, cvtop2(stoi(s),x)));
3081 : }
3082 :
3083 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3084 28 : gel(z,1) = ZX_copy(gel(x,1));
3085 28 : gel(z,2) = gdivgs(gel(x,2),s);
3086 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3087 :
3088 37282 : case t_POLMOD:
3089 37282 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3090 :
3091 91 : case t_RFRAC:
3092 91 : if (s == 1) return gcopy(x);
3093 84 : else if (s == -1) return gneg(x);
3094 84 : return div_rfrac_scal(x, stoi(s));
3095 :
3096 156742 : case t_POL: pari_APPLY_pol_normalized(gdivgs(gel(x,i),s));
3097 0 : case t_SER: pari_APPLY_ser_normalized(gdivgs(gel(x,i),s));
3098 407621 : case t_VEC:
3099 : case t_COL:
3100 895033 : case t_MAT: pari_APPLY_same(gdivgs(gel(x,i),s));
3101 : }
3102 0 : pari_err_TYPE2("/",x, stoi(s));
3103 : return NULL; /* LCOV_EXCL_LINE */
3104 : }
3105 :
3106 : GEN
3107 60481378 : gdivgu(GEN x, ulong s)
3108 : {
3109 60481378 : long tx = typ(x), i;
3110 : pari_sp av;
3111 : GEN z;
3112 :
3113 60481378 : if (!s)
3114 : {
3115 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3116 0 : pari_err_INV("gdivgu",gen_0);
3117 : }
3118 60481352 : switch(tx)
3119 : {
3120 23643111 : case t_INT: return Qdiviu(x, s);
3121 12901775 : case t_REAL: return divru(x,s);
3122 :
3123 210315 : case t_INTMOD:
3124 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3125 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3126 :
3127 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3128 :
3129 727843 : case t_FRAC: z = cgetg(3, t_FRAC);
3130 727843 : i = ugcdiu(gel(x,1), s);
3131 727843 : if (i == 1)
3132 : {
3133 534216 : gel(z,2) = mului(s, gel(x,2));
3134 534216 : gel(z,1) = icopy(gel(x,1));
3135 : }
3136 : else
3137 : {
3138 193627 : gel(z,2) = mului(s/i, gel(x,2));
3139 193627 : gel(z,1) = divis(gel(x,1), i);
3140 : }
3141 727843 : normalize_frac(z);
3142 727843 : fix_frac_if_int(z); return z;
3143 :
3144 11347186 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3145 11347186 : gel(z,1) = gdivgu(gel(x,1),s);
3146 11347186 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3147 :
3148 11551638 : case t_PADIC: /* divpT */
3149 : {
3150 11551638 : GEN p = padic_p(x);
3151 11551638 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3152 11416039 : av = avma;
3153 11416039 : return gc_upto(av, divpp(x, cvtop2(utoi(s),x)));
3154 : }
3155 :
3156 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3157 0 : gel(z,1) = ZX_copy(gel(x,1));
3158 0 : gel(z,2) = gdivgu(gel(x,2),s);
3159 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3160 :
3161 1456 : case t_POLMOD:
3162 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3163 :
3164 56 : case t_RFRAC:
3165 56 : if (s == 1) return gcopy(x);
3166 56 : return div_rfrac_scal(x, utoi(s));
3167 :
3168 451122 : case t_POL: pari_APPLY_pol_normalized(gdivgu(gel(x,i),s));
3169 18011 : case t_SER: pari_APPLY_ser_normalized(gdivgu(gel(x,i),s));
3170 336 : case t_VEC:
3171 : case t_COL:
3172 1162 : case t_MAT: pari_APPLY_same(gdivgu(gel(x,i),s));
3173 : }
3174 0 : pari_err_TYPE2("/",x, utoi(s));
3175 : return NULL; /* LCOV_EXCL_LINE */
3176 : }
3177 :
3178 : /* x / (i*(i+1)) */
3179 : GEN
3180 213819991 : divrunextu(GEN x, ulong i)
3181 : {
3182 213819991 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3183 0 : return divri(x, muluu(i , i+1));
3184 : else
3185 213819991 : return divru(x, i*(i+1));
3186 : }
3187 : /* x / (i*(i+1)) */
3188 : GEN
3189 1583803 : gdivgunextu(GEN x, ulong i)
3190 : {
3191 1583803 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3192 0 : return gdivgu(x, i*(i+1));
3193 : else
3194 1583803 : return gdiv(x, muluu(i, i+1));
3195 : }
3196 :
3197 : /* True shift (exact multiplication by 2^n) */
3198 : GEN
3199 143710503 : gmul2n(GEN x, long n)
3200 : {
3201 : GEN z, a, b;
3202 : long k, l;
3203 :
3204 143710503 : switch(typ(x))
3205 : {
3206 43042993 : case t_INT:
3207 43042993 : if (n>=0) return shifti(x,n);
3208 7637460 : if (!signe(x)) return gen_0;
3209 5082936 : l = vali(x); n = -n;
3210 5082958 : if (n<=l) return shifti(x,-n);
3211 350338 : z = cgetg(3,t_FRAC);
3212 350338 : gel(z,1) = shifti(x,-l);
3213 350338 : gel(z,2) = int2n(n-l); return z;
3214 :
3215 64408719 : case t_REAL:
3216 64408719 : return shiftr(x,n);
3217 :
3218 181078 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3219 181078 : z = cgetg(3,t_INTMOD);
3220 181078 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3221 82917 : gel(z,2) = gc_INT((pari_sp)z, modii(shifti(a,n), b));
3222 82917 : gel(z,1) = icopy(b); return z;
3223 :
3224 217420 : case t_FFELT: return FF_mul2n(x,n);
3225 :
3226 1128845 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3227 1128845 : l = vali(a);
3228 1128845 : k = vali(b);
3229 1128845 : if (n+l >= k)
3230 : {
3231 323410 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3232 229543 : l = n-k; k = -k;
3233 : }
3234 : else
3235 : {
3236 805435 : k = -(l+n); l = -l;
3237 : }
3238 1034978 : z = cgetg(3,t_FRAC);
3239 1034978 : gel(z,1) = shifti(a,l);
3240 1034978 : gel(z,2) = shifti(b,k); return z;
3241 :
3242 11102396 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3243 11102388 : gel(z,1) = gmul2n(gel(x,1),n);
3244 11102384 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3245 :
3246 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3247 105 : gel(z,1) = ZX_copy(gel(x,1));
3248 105 : gel(z,2) = gmul2n(gel(x,2),n);
3249 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3250 :
3251 177009 : case t_POLMOD:
3252 177009 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3253 :
3254 6752292 : case t_POL:
3255 24031981 : pari_APPLY_pol(gmul2n(gel(x,i),n));
3256 103491 : case t_SER:
3257 103491 : if (ser_isexactzero(x)) return gcopy(x);
3258 636895 : pari_APPLY_ser(gmul2n(gel(x,i),n));
3259 16589949 : case t_VEC: case t_COL: case t_MAT:
3260 61024654 : pari_APPLY_same(gmul2n(gel(x,i),n));
3261 :
3262 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3263 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3264 :
3265 7644 : case t_PADIC: /* int2n wrong if n < 0 */
3266 7644 : return gmul(gmul2n(gen_1,n),x);
3267 : }
3268 0 : pari_err_TYPE("gmul2n",x);
3269 : return NULL; /* LCOV_EXCL_LINE */
3270 : }
3271 :
3272 : /*******************************************************************/
3273 : /* */
3274 : /* INVERSE */
3275 : /* */
3276 : /*******************************************************************/
3277 : static GEN
3278 189697 : inv_polmod(GEN T, GEN x)
3279 : {
3280 189697 : GEN z = cgetg(3,t_POLMOD), a;
3281 189697 : gel(z,1) = RgX_copy(T);
3282 189697 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3283 64987 : a = ginv(x);
3284 : else
3285 : {
3286 124710 : if (lg(T) == 5) /* quadratic fields */
3287 13125 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3288 : else
3289 111585 : a = RgXQ_inv(x, T);
3290 : }
3291 189695 : gel(z,2) = a; return z;
3292 : }
3293 :
3294 : GEN
3295 36545700 : ginv(GEN x)
3296 : {
3297 : long s;
3298 : pari_sp av;
3299 : GEN z, y;
3300 :
3301 36545700 : switch(typ(x))
3302 : {
3303 9457058 : case t_INT:
3304 9457058 : if (is_pm1(x)) return icopy(x);
3305 7907597 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3306 7907583 : z = cgetg(3,t_FRAC);
3307 7912170 : gel(z,1) = s<0? gen_m1: gen_1;
3308 7912170 : gel(z,2) = absi(x); return z;
3309 :
3310 10881281 : case t_REAL: return invr(x);
3311 :
3312 6419 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3313 6419 : gel(z,1) = icopy(gel(x,1));
3314 6419 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3315 :
3316 532832 : case t_FRAC: {
3317 532832 : GEN a = gel(x,1), b = gel(x,2);
3318 532832 : s = signe(a);
3319 532832 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3320 247324 : z = cgetg(3,t_FRAC);
3321 247324 : gel(z,1) = icopy(b);
3322 247324 : gel(z,2) = icopy(a);
3323 247324 : normalize_frac(z); return z;
3324 : }
3325 9617818 : case t_COMPLEX:
3326 9617818 : av = avma;
3327 9617818 : return gc_upto(av, divcR(conj_i(x), cxnorm(x)));
3328 :
3329 273 : case t_QUAD:
3330 273 : av = avma;
3331 273 : return gc_upto(av, gdiv(conj_i(x), quadnorm(x)));
3332 :
3333 358358 : case t_PADIC:
3334 : {
3335 358358 : GEN p = padic_p(x), pd = padic_pd(x), u = padic_u(x);
3336 358358 : long d = precp(x);
3337 358358 : if (!signe(u)) pari_err_INV("ginv",x);
3338 358351 : retmkpadic(Zp_inv(u, p, d), icopy(p), icopy(pd), -valp(x), d);
3339 : }
3340 :
3341 189697 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3342 14056 : case t_FFELT: return FF_inv(x);
3343 5322058 : case t_POL: return gred_rfrac_simple(gen_1,x);
3344 26985 : case t_SER: return ser_inv(x);
3345 2933 : case t_RFRAC:
3346 : {
3347 2933 : GEN n = gel(x,1), d = gel(x,2);
3348 2933 : pari_sp av = avma, ltop;
3349 2933 : if (gequal0(n)) pari_err_INV("ginv",x);
3350 :
3351 2933 : n = simplify_shallow(n);
3352 2933 : if (typ(n) != t_POL || varn(n) != varn(d))
3353 : {
3354 2933 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3355 679 : ltop = avma;
3356 679 : z = RgX_Rg_div(d,n);
3357 : } else {
3358 0 : ltop = avma;
3359 0 : z = cgetg(3,t_RFRAC);
3360 0 : gel(z,1) = RgX_copy(d);
3361 0 : gel(z,2) = RgX_copy(n);
3362 : }
3363 679 : stackdummy(av, ltop);
3364 679 : return z;
3365 : }
3366 :
3367 21 : case t_VEC: if (!is_ext_qfr(x)) break;
3368 : case t_QFB:
3369 98200 : return qfbpow(x, gen_m1);
3370 38932 : case t_MAT:
3371 38932 : y = RgM_inv(x);
3372 38925 : if (!y) pari_err_INV("ginv",x);
3373 38855 : return y;
3374 28 : case t_VECSMALL:
3375 : {
3376 28 : long i, lx = lg(x)-1;
3377 28 : y = zero_zv(lx);
3378 112 : for (i=1; i<=lx; i++)
3379 : {
3380 84 : long xi = x[i];
3381 84 : if (xi<1 || xi>lx || y[xi])
3382 0 : pari_err_TYPE("ginv [not a permutation]", x);
3383 84 : y[xi] = i;
3384 : }
3385 28 : return y;
3386 : }
3387 : }
3388 6 : pari_err_TYPE("inverse",x);
3389 : return NULL; /* LCOV_EXCL_LINE */
3390 : }
|
