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 105 : warn_coercion(GEN x, GEN y, GEN z)
38 : {
39 105 : if (DEBUGLEVEL)
40 56 : pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
41 105 : }
42 :
43 : static long
44 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 307125365 : addRc(GEN x, GEN y) {
61 307125365 : GEN z = cgetg(3,t_COMPLEX);
62 307100717 : gel(z,1) = gadd(x,gel(y,1));
63 307093226 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 354208300 : mulRc(GEN x, GEN y)
67 : {
68 354208300 : GEN a = gel(y,1), b = gel(y,2), z = cgetg(3,t_COMPLEX);
69 354223129 : gel(z,1) = (typ(x) != t_INTMOD && isintzero(a))? gen_0: gmul(x, a);
70 354695396 : gel(z,2) = gmul(x, b); return z;
71 : }
72 : static GEN
73 2282452 : divRc(GEN x, GEN y)
74 : {
75 2282452 : GEN a = gel(y,1), b = gel(y,2), z = cgetg(3,t_COMPLEX);
76 2282452 : pari_sp av = avma, av2;
77 2282452 : 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 2233157 : GEN t = gdiv(x, cxnorm(y)), mb = gneg(b);
85 2233150 : av2 = avma;
86 2233150 : gel(z,1) = gmul(t, a);
87 2233139 : gel(z,2) = gmul(t, mb);
88 2233140 : gc_slice_unsafe(av, av2, z+1, 2);
89 : }
90 2282445 : 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 23352739 : divcR(GEN x, GEN y) {
100 23352739 : GEN z = cgetg(3,t_COMPLEX);
101 23352661 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
102 23352237 : 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 37789988 : mulrfrac(GEN x, GEN y)
126 : {
127 : pari_sp av;
128 37789988 : GEN z, a = gel(y,1), b = gel(y,2);
129 37789988 : if (is_pm1(a)) /* frequent special case */
130 : {
131 12771000 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
132 12769680 : return z;
133 : }
134 25018782 : av = avma;
135 25018782 : 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 4984902 : mulTp(GEN x, GEN y) { pari_sp av = avma;
163 4984902 : 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 1204917 : divpT(GEN x, GEN y) { pari_sp av = avma;
174 1204917 : 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 3344833 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 3344833 : if (lgefint(X) == 3) {
182 3021725 : ulong u = Fl_add(itou(x),itou(y), X[2]);
183 3021725 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
184 : }
185 : else {
186 323108 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
187 323108 : gel(z,2) = gc_INT((pari_sp)z, u);
188 : }
189 3344833 : gel(z,1) = icopy(X); return z;
190 : }
191 : static GEN
192 1027929 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
193 1027929 : if (lgefint(X) == 3) {
194 653895 : ulong u = Fl_sub(itou(x),itou(y), X[2]);
195 653895 : 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 1027929 : gel(z,1) = icopy(X); return z;
202 : }
203 : /* cf add_intmod_same */
204 : static GEN
205 2629099 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
206 2629099 : if (lgefint(X) == 3) {
207 2487915 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
208 2487915 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 141184 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x,y), X) );
212 2629095 : gel(z,1) = icopy(X); return z;
213 : }
214 : /* cf add_intmod_same */
215 : static GEN
216 341913 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
217 : {
218 341913 : if (lgefint(X) == 3) {
219 319389 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
220 319326 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
221 : }
222 : else
223 22524 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
224 341851 : 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 9959379 : rfrac_denom_mul_scal(GEN d, GEN y)
237 : {
238 9959379 : GEN D = RgX_Rg_mul(d, y);
239 9959379 : 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 9959379 : return D;
245 : }
246 :
247 : static int
248 57792547 : Leading_is_neg(GEN x)
249 : {
250 122218913 : while (typ(x) == t_POL) x = leading_coeff(x);
251 57792547 : return is_real_t(typ(x))? gsigne(x) < 0: 0;
252 : }
253 :
254 : static int
255 154401681 : 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 57796615 : gred_rfrac_simple(GEN n, GEN d)
260 : {
261 : GEN _1n, _1d, c, cn, cd, z;
262 57796615 : long dd = degpol(d);
263 :
264 57796615 : if (dd <= 0)
265 : {
266 4068 : if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
267 4068 : n = gdiv(n, gel(d,2));
268 4068 : if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
269 4068 : return n;
270 : }
271 57792547 : if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
272 57792547 : _1n = Rg_get_1(n);
273 57792547 : _1d = Rg_get_1(d);
274 57792547 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
275 57792547 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
276 57792540 : cd = content(d);
277 59662490 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
278 57792540 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
279 57792540 : 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 51205659 : if (!gequal1(cn))
296 : {
297 3341395 : if (gequal0(cn)) {
298 987 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
299 21 : c = gen_1;
300 : } else {
301 3340408 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
302 3340408 : c = cn;
303 : }
304 : } else {
305 47864264 : GEN y = cgetg(3,t_RFRAC);
306 47864264 : gel(y,1) = gcopy(n);
307 47864264 : gel(y,2) = RgX_copy(d); return y;
308 : }
309 : }
310 :
311 9927240 : 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 9014624 : cn = numer_i(c);
321 9014624 : cd = denom_i(c);
322 : }
323 9927240 : z = cgetg(3,t_RFRAC);
324 9927240 : gel(z,1) = gmul(n, cn);
325 9927240 : 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 9927240 : if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
328 9927240 : 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 207739 : fix_rfrac(GEN x, long d)
334 : {
335 : GEN z, N, D;
336 207739 : if (!d || typ(x) == t_POL) return x;
337 165879 : z = cgetg(3, t_RFRAC);
338 165879 : N = gel(x,1);
339 165879 : D = gel(x,2);
340 165879 : if (d > 0) {
341 6272 : gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
342 160454 : : monomialcopy(N,d,varn(D));
343 160391 : gel(z, 2) = RgX_copy(D);
344 : } else {
345 5488 : gel(z, 1) = gcopy(N);
346 5488 : gel(z, 2) = RgX_shift(D, -d);
347 : }
348 165879 : return z;
349 : }
350 :
351 : /* assume d != 0 */
352 : static GEN
353 44825494 : gred_rfrac2(GEN n, GEN d)
354 : {
355 : GEN y, z, _1n, _1d;
356 : long v, vd, vn;
357 :
358 44825494 : n = simplify_shallow(n);
359 44825494 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
360 37968143 : d = simplify_shallow(d);
361 37968143 : if (typ(d) != t_POL) return gdiv(n,d);
362 36791247 : vd = varn(d);
363 36791247 : if (typ(n) != t_POL)
364 : {
365 20656203 : if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
366 20654796 : if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
367 0 : pari_err_BUG("gred_rfrac2 [incompatible variables]");
368 : }
369 16135044 : vn = varn(n);
370 16135044 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
371 15993093 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
372 15833320 : _1n = Rg_get_1(n);
373 15833320 : _1d = Rg_get_1(d);
374 15833320 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
375 15833320 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
376 :
377 : /* now n and d are t_POLs in the same variable */
378 15833320 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
379 15833320 : if (!degpol(d))
380 : {
381 12797935 : n = RgX_Rg_div(n,gel(d,2));
382 12797935 : return v? RgX_mulXn(n,v): n;
383 : }
384 :
385 : /* X does not divide gcd(n,d), deg(d) > 0 */
386 3035385 : if (!isinexact(n) && !isinexact(d))
387 : {
388 3035161 : y = RgX_divrem(n, d, &z);
389 3035161 : if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
390 207515 : z = RgX_gcd(d, z);
391 207515 : if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
392 : }
393 207739 : 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 141916672 : Qdivii(GEN x, GEN y)
399 : {
400 141916672 : pari_sp av = avma;
401 : GEN r, q;
402 :
403 141916672 : if (lgefint(y) == 3)
404 : {
405 123562686 : q = Qdiviu(x, y[2]);
406 123564963 : if (signe(y) > 0) return q;
407 11096388 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
408 11096468 : return q;
409 : }
410 18353986 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
411 18355346 : if (equali1(x))
412 : {
413 5186683 : if (!signe(y)) pari_err_INV("gdiv",y);
414 5186578 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
415 : }
416 13168684 : q = dvmdii(x,y,&r);
417 13168537 : if (r == gen_0) return q; /* gen_0 intended */
418 8734368 : r = gcdii(y, r);
419 8734396 : if (lgefint(r) == 3)
420 : {
421 7829711 : ulong t = r[2];
422 7829711 : set_avma(av);
423 7829679 : if (t == 1) q = mkfraccopy(x,y);
424 : else
425 : {
426 2830097 : q = cgetg(3,t_FRAC);
427 2830252 : gel(q,1) = diviuexact(x,t);
428 2830150 : gel(q,2) = diviuexact(y,t);
429 : }
430 : }
431 : else
432 : { /* rare: r and q left on stack for efficiency */
433 904685 : q = cgetg(3,t_FRAC);
434 904699 : gel(q,1) = diviiexact(x,r);
435 904698 : gel(q,2) = diviiexact(y,r);
436 : }
437 8734299 : normalize_frac(q); return q;
438 : }
439 :
440 : static GEN
441 64527959 : 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 148676809 : Qdiviu(GEN x, ulong y)
446 : {
447 : pari_sp av;
448 : ulong r, t;
449 : long s;
450 : GEN q;
451 :
452 148676809 : if (y == 1) return icopy(x);
453 126774802 : if (!y) pari_err_INV("gdiv",gen_0);
454 126778843 : s = signe(x); if (!s) return gen_0;
455 80879979 : if (lgefint(x) == 3)
456 : {
457 75054193 : ulong xx = uel(x,2), g;
458 75054193 : if (xx == 1) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
459 64526923 : g = ugcd(xx, y);
460 64528396 : if (g != 1)
461 : {
462 43169884 : xx /= g; if (y == g) return utoi_sign(xx, s);
463 16959408 : y /= g;
464 : }
465 38317920 : retmkfrac(utoi_sign(xx, s), utoipos(y));
466 : }
467 5825786 : av = avma; q = absdiviu_rem(x,y,&r);
468 5826166 : if (!r)
469 : {
470 2006588 : if (s < 0) togglesign(q);
471 2006588 : return q;
472 : }
473 3819578 : t = ugcd(y, r); set_avma(av);
474 3819587 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
475 1426978 : retmkfrac(diviuexact(x,t), utoipos(y / t));
476 : }
477 :
478 : /* x t_INT, return x/y in reduced form */
479 : GEN
480 1592020 : Qdivis(GEN x, long y)
481 : {
482 : GEN q;
483 1592020 : if (y >= 0) return Qdiviu(x, y);
484 180852 : q = Qdiviu(x, -y);
485 180852 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
486 180852 : 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 31489347 : conj_i(GEN x)
522 : {
523 31489347 : switch(typ(x))
524 : {
525 7459318 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
526 7459318 : return x;
527 :
528 23865910 : 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 31647 : case t_SER: pari_APPLY_ser_normalized(conj_i(gel(x,i)));
540 :
541 153443 : case t_RFRAC:
542 : case t_VEC:
543 : case t_COL:
544 553554 : 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 24592887 : mkpadic_mod(GEN u, GEN p, GEN pd, long e, long d)
646 24592887 : { retmkpadic(modii(u, pd), icopy(p), icopy(pd), e, d); }
647 :
648 : /* x, y compatible PADIC, op = add or sub */
649 : static GEN
650 19971993 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
651 : {
652 19971993 : pari_sp av = avma;
653 : long d, e, r, rx, ry;
654 19971993 : GEN u, z, mod, pdx, pdy, ux, uy, p = padic_p(x);
655 : int swap;
656 :
657 19971993 : e = valp(x);
658 19971993 : r = valp(y); d = r-e;
659 19971993 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
660 19971993 : pdx = padic_pd(x); ux = padic_u(x);
661 19971993 : pdy = padic_pd(y); uy = padic_u(y);
662 19971993 : rx = precp(x);
663 19971993 : ry = precp(y);
664 19971993 : if (d) /* v(x) < v(y) */
665 : {
666 10703644 : r = d+ry; z = powiu(p,d);
667 10703711 : if (r < rx) mod = mulii(z,pdy); else { r = rx; mod = pdx; }
668 10703738 : z = mulii(z, uy);
669 10703711 : u = swap? op(z, ux): op(ux, z);
670 : }
671 : else
672 : {
673 : long c;
674 9268349 : if (ry < rx) { r=ry; mod = pdy; } else { r=rx; mod = pdx; }
675 9268349 : u = swap? op(uy, ux): op(ux, uy);
676 9269334 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
677 : {
678 138271 : set_avma(av); return zeropadic(p, e+r);
679 : }
680 9130945 : if (c)
681 : {
682 3351557 : mod = diviiexact(mod, powiu(p,c));
683 3351556 : r -= c;
684 3351556 : e += c;
685 : }
686 : }
687 19834602 : 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 4762937 : addQp(GEN x, GEN y)
696 : {
697 4762937 : pari_sp av = avma;
698 4762937 : long d, r, e, vy = valp(y), py = precp(y);
699 4762937 : GEN q, mod, u, p = padic_p(y);
700 :
701 4762937 : e = Q_pvalrem(x, p, &x);
702 4762930 : d = vy - e; r = d + py;
703 4762930 : if (r <= 0) { set_avma(av); return gcopy(y); }
704 4761044 : mod = padic_pd(y);
705 4761044 : u = padic_u(y);
706 4761044 : if (d > 0)
707 : {
708 1376545 : q = powiu(p,d);
709 1376549 : mod = mulii(mod, q);
710 1376553 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
711 1376553 : u = addii(x, mulii(u, q));
712 : }
713 3384499 : 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 2979167 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
724 2979167 : u = addii(u, x);
725 2979160 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
726 : {
727 1205 : set_avma(av); return zeropadic(p,e+r);
728 : }
729 2977954 : if (c)
730 : {
731 970644 : mod = diviiexact(mod, powiu(p,c));
732 970643 : r -= c;
733 970643 : e += c;
734 : }
735 : }
736 4759843 : 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 13455135 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
763 13455135 : { 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 408206 : add_ser_scal(GEN y, GEN x)
768 : {
769 : long i, v, ly, vy;
770 : GEN z;
771 :
772 408206 : if (isrationalzero(x)) return gcopy(y);
773 381921 : ly = lg(y);
774 381921 : v = valser(y);
775 381921 : if (v < 3-ly) return gcopy(y);
776 : /* v + ly >= 3 */
777 381669 : if (v < 0)
778 : {
779 1162 : z = cgetg(ly,t_SER); z[1] = y[1];
780 3276 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
781 1162 : gel(z,i) = gadd(x,gel(y,i)); i++;
782 3157 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
783 1162 : return normalizeser(z);
784 : }
785 380507 : vy = varn(y);
786 380507 : if (v > 0)
787 : {
788 19355 : if (ser_isexactzero(y))
789 9394 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
790 9961 : y -= v; ly += v;
791 9961 : z = cgetg(ly,t_SER);
792 9961 : x = gcopy(x);
793 20559 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
794 : }
795 361152 : else if (ser_isexactzero(y)) return gcopy(y);
796 : else
797 : { /* v = 0, ly >= 3 */
798 361145 : z = cgetg(ly,t_SER);
799 361145 : x = gadd(x, gel(y,2));
800 361145 : i = 3;
801 : }
802 1599007 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
803 371106 : gel(z,2) = x;
804 371106 : z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
805 371106 : return gequal0(x)? normalizeser(z): z;
806 : }
807 : static long
808 7162062 : _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 3581430 : ser_add(GEN x, GEN y)
812 : {
813 3581430 : long i, lx,ly, n = valser(y) - valser(x);
814 : GEN z;
815 3581430 : if (n < 0) { n = -n; swap(x,y); }
816 : /* valser(x) <= valser(y) */
817 3581430 : lx = _serprec(x);
818 3581430 : 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 3580632 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
824 3580632 : if (n)
825 : {
826 107803 : if (n+2 > lx) return gcopy(x);
827 107103 : z = cgetg(ly,t_SER);
828 819228 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
829 506771 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
830 : } else {
831 3472829 : z = cgetg(ly,t_SER);
832 17589564 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
833 : }
834 3579932 : z[1] = x[1]; return normalizeser(z);
835 : }
836 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
837 : static GEN
838 8811954 : add_rfrac_scal(GEN y, GEN x)
839 : {
840 : pari_sp av;
841 : GEN n;
842 :
843 8811954 : if (isintzero(x)) return gcopy(y); /* frequent special case */
844 5138189 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
845 5138189 : 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 42804949 : add_scal(GEN y, GEN x, long ty)
851 : {
852 42804949 : switch(ty)
853 : {
854 37882517 : case t_POL: return RgX_Rg_add(y, x);
855 408178 : case t_SER: return add_ser_scal(y, x);
856 4468855 : case t_RFRAC: return add_rfrac_scal(y, x);
857 0 : case t_COL: return RgC_Rg_add(y, x);
858 45387 : case t_VEC:
859 45387 : 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 14998135 : setfrac(GEN z, GEN a, GEN b)
869 : {
870 14998135 : gel(z,1) = icopy_avma(a, (pari_sp)z);
871 14998130 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
872 14998176 : set_avma((pari_sp)gel(z,2)); return z;
873 : }
874 : /* z <- a / (b*Q), (Q,a) = 1 */
875 : static GEN
876 14120186 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
877 : {
878 14120186 : GEN q = Qdivii(a, b); /* != 0 */
879 14120240 : if (typ(q) == t_INT)
880 : {
881 1904621 : gel(z,1) = gc_INT((pari_sp)Q, q);
882 1904621 : gel(z,2) = Q; return z;
883 : }
884 12215619 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
885 : }
886 : static GEN
887 40597181 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
888 : {
889 40597181 : GEN x1 = gel(x,1), x2 = gel(x,2);
890 40597181 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
891 40597181 : int s = cmpii(x2, y2);
892 :
893 : /* common denominator: (x1 op y1) / x2 */
894 40597183 : if (!s)
895 : {
896 15298988 : pari_sp av = avma;
897 15298988 : return gc_upto(av, Qdivii(op(x1, y1), x2));
898 : }
899 25298195 : z = cgetg(3, t_FRAC);
900 25298205 : if (s < 0)
901 : {
902 17685313 : Q = dvmdii(y2, x2, &r);
903 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
904 17685270 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
905 9723011 : delta = gcdii(x2,r);
906 : }
907 : else
908 : {
909 7612892 : Q = dvmdii(x2, y2, &r);
910 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
911 7612908 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
912 1454934 : delta = gcdii(y2,r);
913 : }
914 : /* delta = gcd(x2,y2) */
915 11177978 : if (equali1(delta))
916 : { /* numerator is nonzero */
917 8395412 : gel(z,1) = gc_INT((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
918 8395412 : gel(z,2) = mulii(x2,y2); return z;
919 : }
920 2782559 : x2 = diviiexact(x2,delta);
921 2782564 : y2 = diviiexact(y2,delta);
922 2782564 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
923 2782567 : q = dvmdii(n, delta, &r);
924 2782565 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
925 2528582 : r = gcdii(delta, r);
926 2528587 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
927 2528587 : 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 3018271 : add_rfrac(GEN x, GEN y)
933 : {
934 3018271 : pari_sp av = avma;
935 3018271 : GEN x1 = gel(x,1), x2 = gel(x,2);
936 3018271 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
937 :
938 3018271 : delta = RgX_gcd(x2,y2);
939 3018271 : 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 3017627 : x2 = RgX_div(x2,delta);
946 3017627 : y2 = RgX_div(y2,delta);
947 3017627 : n = gadd(gmul(x1,y2), gmul(y1,x2));
948 3017627 : 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 2296209 : if (degpol(n) == 0)
959 1150575 : return gc_upto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
960 1145634 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
961 1145634 : if (isexactzero(r))
962 : {
963 : GEN z;
964 227982 : d = RgX_mul(x2, y2);
965 : /* "constant" denominator ? */
966 227982 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
967 227982 : return gc_upto(av, z);
968 : }
969 917652 : r = RgX_gcd(delta, r);
970 917652 : if (degpol(r))
971 : {
972 160642 : n = RgX_div(n, r);
973 160642 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
974 : }
975 : else
976 757010 : d = RgX_mul(gel(x,2), y2);
977 917652 : return gc_upto(av, gred_rfrac_simple(n, d));
978 : }
979 :
980 : GEN
981 6138460880 : gadd(GEN x, GEN y)
982 : {
983 6138460880 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
984 : pari_sp av;
985 : GEN z, p1;
986 :
987 6138460880 : if (tx == ty) switch(tx) /* shortcut to generic case */
988 : {
989 3055803849 : case t_INT: return addii(x,y);
990 1657339118 : case t_REAL: return addrr(x,y);
991 1176494 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
992 1176494 : z = cgetg(3,t_INTMOD);
993 1176494 : if (X==Y || equalii(X,Y))
994 1176480 : 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 34504400 : case t_FRAC: return addsub_frac(x,y,addii);
1001 360011175 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1002 359526876 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1003 360287108 : 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 359583355 : gel(z,1) = gadd(gel(x,1),gel(y,1));
1009 359706826 : return z;
1010 14529620 : case t_PADIC:
1011 14529620 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1012 14529387 : 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 10151518 : case t_POLMOD:
1019 10151518 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1020 10144350 : 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 13133927 : case t_FFELT: return FF_add(x,y);
1023 138941130 : case t_POL:
1024 138941130 : vx = varn(x);
1025 138941130 : vy = varn(y);
1026 138941130 : if (vx != vy) {
1027 4863955 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1028 1801522 : else return RgX_Rg_add(y, x);
1029 : }
1030 134077175 : return RgX_add(x, y);
1031 3576068 : case t_SER:
1032 3576068 : vx = varn(x);
1033 3576068 : vy = varn(y);
1034 3576068 : 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 3576040 : return ser_add(x, y);
1039 4310172 : case t_RFRAC:
1040 4310172 : vx = varn(gel(x,2));
1041 4310172 : vy = varn(gel(y,2));
1042 4310172 : 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 3018271 : return add_rfrac(x,y);
1047 2115703 : case t_VEC:
1048 2115703 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1049 2115703 : return RgV_add(x,y);
1050 1113429 : case t_COL:
1051 1113429 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1052 1113429 : 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 845240248 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1064 :
1065 845240248 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1066 : {
1067 730350371 : case t_INT:
1068 : switch(ty)
1069 : {
1070 441113703 : case t_REAL: return addir(x,y);
1071 2138731 : case t_INTMOD:
1072 2138731 : z = cgetg(3, t_INTMOD);
1073 2138731 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1074 27563857 : case t_FRAC: z = cgetg(3,t_FRAC);
1075 27563853 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1076 27563791 : gel(z,2) = icopy(gel(y,2)); return z;
1077 252614085 : case t_COMPLEX: return addRc(x, y);
1078 5597450 : case t_PADIC:
1079 5597450 : if (!signe(x)) return gcopy(y);
1080 4533165 : return addQp(x,y);
1081 1113 : case t_QUAD: return addRq(x, y);
1082 1365080 : case t_FFELT: return FF_Z_add(y,x);
1083 : }
1084 :
1085 : case t_REAL:
1086 58604471 : switch(ty)
1087 : {
1088 14402205 : case t_FRAC:
1089 14402205 : if (!signe(gel(y,1))) return rcopy(x);
1090 14402205 : 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 14390125 : av = avma; z = addir(gel(y,1), mulir(gel(y,2),x));
1096 14388583 : return gc_leaf(av, divri(z,gel(y,2)));
1097 44202196 : 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 10327465 : 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 0 : 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 48496502 : switch(ty)
1152 : {
1153 7693 : case t_MAT:
1154 7693 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1155 7693 : if (isrationalzero(x)) return gcopy(y);
1156 7616 : return RgM_Rg_add(y, x);
1157 186395 : case t_COL:
1158 186395 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1159 186395 : return RgC_Rg_add(y, x);
1160 2409752 : case t_POLMOD: /* is_const_t(tx) in this case */
1161 2409752 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1162 : }
1163 45892662 : if (is_scalar_t(tx)) {
1164 42812787 : if (tx == t_POLMOD)
1165 : {
1166 117809 : vx = varn(gel(x,1));
1167 117809 : vy = gvar(y);
1168 117809 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1169 : else
1170 87758 : 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 42694978 : return add_scal(y, x, ty);
1174 : }
1175 : /* x and y are not scalars, ty != t_MAT */
1176 3079867 : vx = gvar(x);
1177 3079872 : vy = gvar(y);
1178 3079872 : 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 3057113 : switch(tx)
1185 : {
1186 3056616 : 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 3051198 : 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 270736848 : gaddsg(long x, GEN y)
1223 : {
1224 270736848 : long ty = typ(y);
1225 : GEN z;
1226 :
1227 270736848 : switch(ty)
1228 : {
1229 124844351 : case t_INT: return addsi(x,y);
1230 119268577 : 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 15248685 : case t_FRAC: z = cgetg(3,t_FRAC);
1235 15248685 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1236 15248685 : gel(z,2) = icopy(gel(y,2)); return z;
1237 8225330 : case t_COMPLEX:
1238 8225330 : z = cgetg(3, t_COMPLEX);
1239 8225330 : gel(z,1) = gaddsg(x, gel(y,1));
1240 8225330 : gel(z,2) = gcopy(gel(y,2)); return z;
1241 :
1242 3137846 : default: return gadd(stoi(x), y);
1243 : }
1244 : }
1245 :
1246 : GEN
1247 3238198 : gsubsg(long x, GEN y)
1248 : {
1249 : GEN z, a, b;
1250 : pari_sp av;
1251 :
1252 3238198 : switch(typ(y))
1253 : {
1254 276785 : case t_INT: return subsi(x,y);
1255 1234943 : 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 732508 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1260 732508 : gel(z,1) = gc_INT((pari_sp)z, subii(mulis(b,x), a));
1261 732508 : gel(z,2) = icopy(gel(y,2)); return z;
1262 956792 : case t_COMPLEX:
1263 956792 : z = cgetg(3, t_COMPLEX);
1264 956792 : gel(z,1) = gsubsg(x, gel(y,1));
1265 956792 : gel(z,2) = gneg(gel(y,2)); return z;
1266 : }
1267 37114 : av = avma;
1268 37114 : return gc_upto(av, gadd(stoi(x), gneg_i(y)));
1269 : }
1270 :
1271 : /********************************************************************/
1272 : /** **/
1273 : /** SUBTRACTION **/
1274 : /** **/
1275 : /********************************************************************/
1276 :
1277 : GEN
1278 2885304372 : gsub(GEN x, GEN y)
1279 : {
1280 2885304372 : long tx = typ(x), ty = typ(y);
1281 : pari_sp av;
1282 : GEN z;
1283 2885304372 : if (tx == ty) switch(tx) /* shortcut to generic case */
1284 : {
1285 2093825505 : case t_INT: return subii(x,y);
1286 586046602 : case t_REAL: return subrr(x,y);
1287 1027943 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1288 1027943 : z = cgetg(3,t_INTMOD);
1289 1027943 : if (X==Y || equalii(X,Y))
1290 1027929 : 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 6092789 : case t_FRAC: return addsub_frac(x,y, subii);
1297 102429140 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1298 102464652 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1299 102374822 : 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 102356273 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1305 102341449 : return z;
1306 5442665 : case t_PADIC:
1307 5442665 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1308 5442665 : 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 870478 : case t_POLMOD:
1315 870478 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1316 870443 : 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 240357 : case t_FFELT: return FF_sub(x,y);
1319 20224606 : case t_POL: {
1320 20224606 : long vx = varn(x);
1321 20224606 : long vy = varn(y);
1322 20224606 : if (vx != vy) {
1323 1150606 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1324 85352 : else return Rg_RgX_sub(x, y);
1325 : }
1326 19074000 : return RgX_sub(x, y);
1327 : }
1328 299446 : case t_VEC:
1329 299446 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1330 299446 : return RgV_sub(x,y);
1331 3607157 : case t_COL:
1332 3607157 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1333 3607157 : return RgC_sub(x,y);
1334 253977 : case t_MAT: {
1335 253977 : long lx = lg(x);
1336 253977 : if (lg(y) != lx) pari_err_OP("+",x,y);
1337 253977 : if (lx == 1) return cgetg(1, t_MAT);
1338 173885 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1339 173884 : return RgM_sub(x,y);
1340 : }
1341 2448387 : case t_RFRAC: case t_SER: break;
1342 :
1343 0 : default: pari_err_TYPE2("+",x,y);
1344 : }
1345 66552365 : av = avma;
1346 66552365 : return gc_upto(av, gadd(x,gneg_i(y)));
1347 : }
1348 :
1349 : /********************************************************************/
1350 : /** **/
1351 : /** MULTIPLICATION **/
1352 : /** **/
1353 : /********************************************************************/
1354 : static GEN
1355 316386 : mul_ser_scal(GEN x, GEN t)
1356 : {
1357 316386 : if (isexactzero(t)) return gmul(Rg_get_0(x), t);
1358 313138 : if (isint1(t)) return gcopy(x);
1359 262514 : 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 3442782 : 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 10454828 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1370 : {
1371 10454828 : pari_sp av = avma;
1372 : GEN z;
1373 :
1374 10454828 : 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 630 : case t_INTMOD: case t_POLMOD:
1382 630 : n = gmul(n, x);
1383 630 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1384 630 : return gc_upto(av, gdiv(n,d));
1385 : }
1386 10454177 : z = gred_rfrac2(x, d);
1387 10454177 : n = simplify_shallow(n);
1388 10454177 : if (typ(z) == t_RFRAC)
1389 : {
1390 7929477 : n = gmul(gel(z,1), n);
1391 7929477 : d = gel(z,2);
1392 7929477 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1393 0 : z = RgX_Rg_div(n, d);
1394 : else
1395 7929477 : z = gred_rfrac_simple(n, d);
1396 : }
1397 : else
1398 2524700 : z = gmul(z, n);
1399 10454177 : return gc_upto(av, z);
1400 : }
1401 : static GEN
1402 154324227 : mul_scal(GEN y, GEN x, long ty)
1403 : {
1404 154324227 : switch(ty)
1405 : {
1406 145601259 : case t_POL:
1407 145601259 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1408 133535564 : return RgX_Rg_mul(y, x);
1409 308245 : case t_SER: return mul_ser_scal(y, x);
1410 8414709 : 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 7645330 : mul_self_scal(GEN x, GEN y)
1419 640610857 : { 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 7910285 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1432 : {
1433 : GEN z, X, Y;
1434 7910285 : pari_sp av = avma;
1435 :
1436 7910285 : X = gred_rfrac2(x1, y2);
1437 7910285 : Y = gred_rfrac2(y1, x2);
1438 7910285 : if (typ(X) == t_RFRAC)
1439 : {
1440 6631189 : if (typ(Y) == t_RFRAC) {
1441 6565032 : x1 = gel(X,1);
1442 6565032 : x2 = gel(X,2);
1443 6565032 : y1 = gel(Y,1);
1444 6565032 : y2 = gel(Y,2);
1445 6565032 : 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 7910285 : return gc_upto(av, z);
1454 : }
1455 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1456 : static GEN
1457 262344 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1458 : {
1459 262344 : pari_sp av = avma;
1460 262344 : GEN X = gred_rfrac2(x1, y2);
1461 262344 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1462 : {
1463 255764 : x2 = RgX_mul(gel(X,2), x2);
1464 255764 : x1 = gel(X,1);
1465 : }
1466 : else
1467 6580 : x1 = X;
1468 262344 : return gc_upto(av, gred_rfrac_simple(x1, x2));
1469 : }
1470 :
1471 : /* Mod(y, Y) * x, assuming x scalar */
1472 : static GEN
1473 3238383 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1474 3238383 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1475 :
1476 : /* cf mulqq */
1477 : static GEN
1478 5860097 : quad_polmod_mul(GEN T, GEN x, GEN y)
1479 : {
1480 5860097 : GEN b = gel(T,3), c = gel(T,2);
1481 5860097 : 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 5860097 : z = cgetg(4, t_POL); z[1] = x[1]; av = avma;
1485 5860097 : U = gmul(ux, uy);
1486 5860097 : V = gmul(vx, vy); s = gmul(c, V);
1487 5860097 : E = gmul(gadd(ux, vx), gadd(uy, vy));
1488 5860097 : if (typ(b) != t_INT) F = gadd(U, gmul(gaddgs(b, 1), V)); /* = U + (b+1) V */
1489 : else
1490 : { /* minor optimization */
1491 5860076 : 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 5860097 : av2 = avma;
1498 5860097 : gel(z,2) = gsub(U, s);
1499 5860097 : gel(z,3) = gsub(E, F);
1500 5860097 : gc_slice_unsafe(av, av2, z+2, 2);
1501 5860097 : return normalizepol_lg(z,4);
1502 : }
1503 : /* Mod(x,T) * Mod(y,T) */
1504 : static GEN
1505 8528815 : mul_polmod_same(GEN T, GEN x, GEN y)
1506 : {
1507 8528815 : GEN z = cgetg(3,t_POLMOD), a;
1508 8528815 : long v = varn(T), lx = lg(x), ly = lg(y);
1509 8528815 : gel(z,1) = RgX_copy(T);
1510 : /* x * y mod T optimised */
1511 8528815 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1512 7902082 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1513 1545858 : a = gmul(x, y);
1514 : else
1515 : {
1516 6982957 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1517 5854819 : a = quad_polmod_mul(T, x, y);
1518 : else
1519 1128138 : a = RgXQ_mul(x, y, gel(z,1));
1520 : }
1521 8528815 : 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 11647060 : mulcIR(GEN x, GEN y)
1609 : {
1610 11647060 : GEN z = cgetg(3,t_COMPLEX);
1611 11647217 : pari_sp av = avma;
1612 11647217 : gel(z,1) = gc_upto(av, gneg(gmul(y,gel(x,2))));
1613 11647378 : gel(z,2) = gmul(y, gel(x,1));
1614 11647400 : return z;
1615 :
1616 : }
1617 : /* x,y COMPLEX */
1618 : static GEN
1619 278974715 : mulcc(GEN x, GEN y)
1620 : {
1621 278974715 : GEN xr = gel(x,1), xi = gel(x,2);
1622 278974715 : GEN yr = gel(y,1), yi = gel(y,2);
1623 : GEN p1, p2, p3, p4, z;
1624 : pari_sp tetpil, av;
1625 :
1626 278974715 : if (isintzero(xr))
1627 : {
1628 15645257 : if (isintzero(yr)) {
1629 7122429 : av = avma;
1630 7122429 : return gc_upto(av, gneg(gmul(xi,yi)));
1631 : }
1632 8522654 : return mulcIR(y, xi);
1633 : }
1634 263319278 : if (isintzero(yr)) return mulcIR(x, yi);
1635 :
1636 260185098 : 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 260192615 : 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 259633259 : p1 = gmul(xr,yr);
1697 259583167 : p2 = gmul(xi,yi);
1698 259582018 : p3 = gmul(xr,yi);
1699 259588249 : p4 = gmul(xi,yr);
1700 259592519 : tetpil = avma;
1701 259592519 : gel(z,1) = gsub(p1,p2);
1702 259451748 : gel(z,2) = gadd(p3,p4);
1703 259466151 : 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 259859445 : gc_slice_unsafe(av,tetpil, z+1,2); return z;
1712 : }
1713 : /* x,y PADIC */
1714 : static GEN
1715 17503053 : mulpp(GEN x, GEN y) {
1716 17503053 : GEN pd, p = padic_p(x), ux = padic_u(x), uy = padic_u(y);
1717 17503053 : long e = valp(x) + valp(y), dx, dy;
1718 :
1719 17503053 : if (!equalii(p, padic_p(y))) pari_err_OP("*",x,y);
1720 17503010 : if (!signe(ux) || !signe(uy)) return zeropadic(p, e);
1721 16930340 : dx = precp(x); dy = precp(y);
1722 16930340 : if (dx > dy) { pd = padic_pd(y); dx = dy; } else pd = padic_pd(x);
1723 16930340 : 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 15630268 : mulcxI(GEN x)
1748 : {
1749 : GEN z;
1750 15630268 : switch(typ(x))
1751 : {
1752 1007054 : case t_INT:
1753 1007054 : if (!signe(x)) return gen_0;
1754 : case t_REAL: case t_FRAC:
1755 2611208 : return mkcomplex(gen_0, x);
1756 13018829 : case t_COMPLEX:
1757 13018829 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1758 13005389 : z = cgetg(3,t_COMPLEX);
1759 13005759 : gel(z,1) = gneg(gel(x,2));
1760 13006181 : gel(z,2) = gel(x,1); return z;
1761 56 : default:
1762 56 : return gmul(gen_I(), x);
1763 : }
1764 : }
1765 : GEN
1766 2985911 : mulcxmI(GEN x)
1767 : {
1768 : GEN z;
1769 2985911 : 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 2685092 : case t_COMPLEX:
1776 2685092 : if (isintzero(gel(x,1))) return gel(x,2);
1777 2676676 : z = cgetg(3,t_COMPLEX);
1778 2676645 : gel(z,1) = gel(x,2);
1779 2676645 : 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 5564705 : mulcxpowIs(GEN x, long k)
1787 : {
1788 5564705 : switch (k & 3)
1789 : {
1790 1369270 : case 1: return mulcxI(x);
1791 1362281 : case 2: return gneg(x);
1792 1403440 : case 3: return mulcxmI(x);
1793 : }
1794 1429714 : return x;
1795 : }
1796 :
1797 : /* caller will assume l > 2 later */
1798 : static GEN
1799 5706368 : init_ser(long l, long v, long e)
1800 : {
1801 5706368 : GEN z = cgetg(l, t_SER);
1802 5706368 : z[1] = evalvalser(e) | evalvarn(v) | evalsigne(1); return z;
1803 : }
1804 :
1805 : /* fill in coefficients of t_SER z from coeffs of t_POL y */
1806 : static GEN
1807 5693993 : fill_ser(GEN z, GEN y)
1808 : {
1809 5693993 : long i, lx = lg(z), ly = lg(y); /* lx > 2 */
1810 5693993 : if (ly >= lx) {
1811 25570893 : for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
1812 : } else {
1813 335415 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1814 239573 : for ( ; i < lx; i++) gel(z,i) = gen_0;
1815 : }
1816 : /* dangerous case (already normalized), don't use normalizeser */
1817 5693993 : if (ly == 3 && !signe(y)) { setsigne(z, 0); return z; }
1818 5693895 : return normalizeser(z);
1819 : }
1820 : /* assume typ(x) = t_VEC */
1821 : static int
1822 98 : is_ext_qfr(GEN x)
1823 91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
1824 189 : && typ(gel(x,2)) == t_REAL; }
1825 :
1826 : GEN
1827 8688180864 : gmul(GEN x, GEN y)
1828 : {
1829 : long tx, ty, lx, ly, vx, vy, i, l;
1830 : pari_sp av, tetpil;
1831 : GEN z, p1;
1832 :
1833 8688180864 : if (x == y) return gsqr(x);
1834 7749072622 : tx = typ(x); ty = typ(y);
1835 7749072622 : if (tx == ty) switch(tx)
1836 : {
1837 3547377377 : case t_INT: return mulii(x,y);
1838 1894414407 : case t_REAL: return mulrr(x,y);
1839 1286504 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1840 1286504 : z = cgetg(3,t_INTMOD);
1841 1286504 : if (X==Y || equalii(X,Y))
1842 1286469 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1843 35 : gel(z,1) = gcdii(X,Y);
1844 35 : warn_coercion(X,Y,gel(z,1));
1845 35 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1846 35 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
1847 : }
1848 22833166 : case t_FRAC:
1849 : {
1850 22833166 : GEN x1 = gel(x,1), x2 = gel(x,2);
1851 22833166 : GEN y1 = gel(y,1), y2 = gel(y,2);
1852 22833166 : z=cgetg(3,t_FRAC);
1853 22833165 : p1 = gcdii(x1, y2);
1854 22833164 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1855 22833165 : p1 = gcdii(x2, y1);
1856 22833162 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1857 22833162 : tetpil = avma;
1858 22833162 : gel(z,2) = mulii(x2,y2);
1859 22833165 : gel(z,1) = mulii(x1,y1);
1860 22833163 : fix_frac_if_int_GC(z,tetpil); return z;
1861 : }
1862 269730927 : case t_COMPLEX: return mulcc(x, y);
1863 12509565 : case t_PADIC: return mulpp(x, y);
1864 882 : case t_QUAD:
1865 882 : if (!ZX_equal(gel(x,1), gel(y,1))) pari_err_OP("*",x,y);
1866 882 : return mulqq(x, y);
1867 13689074 : case t_FFELT: return FF_mul(x, y);
1868 10950787 : case t_POLMOD:
1869 10950787 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1870 8275842 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1871 2674945 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1872 131559103 : case t_POL:
1873 131559103 : vx = varn(x);
1874 131559103 : vy = varn(y);
1875 131559103 : if (vx != vy) {
1876 46850451 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1877 27724525 : else return RgX_Rg_mul(y, x);
1878 : }
1879 84708652 : return RgX_mul(x, y);
1880 :
1881 4411917 : case t_SER: {
1882 4411917 : vx = varn(x);
1883 4411917 : vy = varn(y);
1884 4411917 : if (vx != vy) {
1885 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1886 3675 : return mul_ser_scal(y, x);
1887 : }
1888 4408242 : lx = minss(lg(x), lg(y));
1889 4408242 : if (lx == 2) return zeroser(vx, valser(x)+valser(y));
1890 4408228 : av = avma; z = init_ser(lx, vx, valser(x)+valser(y));
1891 4408228 : x = ser2pol_i(x, lx);
1892 4408228 : y = ser2pol_i(y, lx);
1893 4408228 : y = RgXn_mul(x, y, lx-2);
1894 4408228 : return gc_GEN(av, fill_ser(z,y));
1895 : }
1896 42 : case t_VEC:
1897 42 : if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
1898 : /* fall through, handle extended t_QFB */
1899 198 : case t_QFB: return qfbcomp(x,y);
1900 6723048 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1901 4126392 : case t_MAT: return RgM_mul(x, y);
1902 :
1903 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1904 1631 : z = cgetg_copy(x, &l);
1905 1631 : if (l != lg(y)) break;
1906 17325 : for (i=1; i<l; i++)
1907 : {
1908 15694 : long yi = y[i];
1909 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1910 15694 : z[i] = x[yi];
1911 : }
1912 1631 : return z;
1913 :
1914 0 : default:
1915 0 : pari_err_TYPE2("*",x,y);
1916 : }
1917 : /* tx != ty */
1918 1831928099 : if (is_const_t(ty) && is_const_t(tx)) {
1919 1668526611 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1920 1668526611 : switch(tx) {
1921 1539071604 : case t_INT:
1922 : switch(ty)
1923 : {
1924 996220809 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1925 1342377 : case t_INTMOD:
1926 1342377 : z = cgetg(3, t_INTMOD);
1927 1342377 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1928 102325633 : case t_FRAC:
1929 102325633 : if (!signe(x)) return gen_0;
1930 85186181 : z=cgetg(3,t_FRAC);
1931 85186381 : p1 = gcdii(x,gel(y,2));
1932 85185434 : if (equali1(p1))
1933 : {
1934 46043912 : set_avma((pari_sp)z);
1935 46043899 : gel(z,2) = icopy(gel(y,2));
1936 46044003 : gel(z,1) = mulii(gel(y,1), x);
1937 : }
1938 : else
1939 : {
1940 39141823 : x = diviiexact(x,p1); tetpil = avma;
1941 39141155 : gel(z,2) = diviiexact(gel(y,2), p1);
1942 39140856 : gel(z,1) = mulii(gel(y,1), x);
1943 39141444 : fix_frac_if_int_GC(z,tetpil);
1944 : }
1945 85186126 : return z;
1946 436284365 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1947 4198903 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1948 1904 : case t_QUAD: return mulRq(x,y);
1949 1607125 : case t_FFELT: return FF_Z_mul(y,x);
1950 : }
1951 :
1952 : case t_REAL:
1953 121831669 : switch(ty)
1954 : {
1955 37789967 : case t_FRAC: return mulrfrac(x, y);
1956 84041660 : case t_COMPLEX: return mulRc(x, y);
1957 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1958 21 : default: pari_err_TYPE2("*",x,y);
1959 : }
1960 :
1961 8036 : case t_INTMOD:
1962 : switch(ty)
1963 : {
1964 7133 : case t_FRAC: { GEN X = gel(x,1);
1965 7133 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1966 7133 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1967 : }
1968 49 : case t_COMPLEX: return mulRc(x,y);
1969 259 : case t_PADIC: { GEN X = gel(x,1);
1970 259 : z = cgetg(3, t_INTMOD);
1971 259 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1972 : }
1973 63 : case t_QUAD: return mulRq(x, y);
1974 532 : case t_FFELT:
1975 532 : if (!equalii(gel(x,1),FF_p_i(y)))
1976 0 : pari_err_OP("*",x,y);
1977 532 : return FF_Z_mul(y,gel(x,2));
1978 : }
1979 :
1980 : case t_FRAC:
1981 : switch(ty)
1982 : {
1983 5372483 : case t_COMPLEX: return mulRc(x, y);
1984 2582744 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
1985 343 : case t_QUAD: return mulRq(x, y);
1986 2079 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
1987 : }
1988 :
1989 : case t_FFELT:
1990 35 : pari_err_TYPE2("*",x,y);
1991 :
1992 21 : case t_COMPLEX:
1993 : switch(ty)
1994 : {
1995 14 : case t_PADIC:
1996 14 : return Zp_nosquare_m1(padic_p(y))? mulRc(y, x): mulTp(x, y);
1997 7 : case t_QUAD:
1998 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
1999 7 : return mulqf(y, x, lx);
2000 : }
2001 :
2002 : case t_PADIC: /* ty == t_QUAD */
2003 28 : return (kro_quad(y,padic_p(x))== -1)? mulRq(x, y): mulTp(y, x);
2004 : }
2005 : }
2006 :
2007 165058443 : if (is_matvec_t(ty)) switch(tx)
2008 : {
2009 122060 : case t_VEC:
2010 : switch(ty) {
2011 14826 : case t_COL: return RgV_RgC_mul(x,y);
2012 107234 : case t_MAT: return RgV_RgM_mul(x,y);
2013 : }
2014 0 : break;
2015 1687 : case t_COL:
2016 : switch(ty) {
2017 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2018 0 : case t_MAT: return RgC_RgM_mul(x,y);
2019 : }
2020 0 : break;
2021 880378 : case t_MAT:
2022 : switch(ty) {
2023 0 : case t_VEC: return RgM_RgV_mul(x,y);
2024 880378 : case t_COL: return RgM_RgC_mul(x,y);
2025 : }
2026 : default:
2027 5519924 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
2028 5519931 : return mul_self_scal(y, x);
2029 : }
2030 161829729 : if (is_matvec_t(tx))
2031 : {
2032 2125170 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2033 2125323 : return mul_self_scal(x, y);
2034 : }
2035 159704817 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2036 : /* tx < ty, !ismatvec(x and y) */
2037 :
2038 159704817 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2039 3215290 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2040 156489527 : if (is_scalar_t(tx)) {
2041 153020309 : if (tx == t_POLMOD) {
2042 3159530 : vx = varn(gel(x,1));
2043 3159530 : vy = gvar(y);
2044 3159530 : if (vx != vy) {
2045 2907313 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2046 23093 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2047 : }
2048 : /* error if ty == t_SER */
2049 252217 : av = avma; y = gmod(y, gel(x,1));
2050 252210 : return gc_upto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2051 : }
2052 149860779 : return mul_scal(y, x, ty);
2053 : }
2054 :
2055 : /* x and y are not scalars, nor matvec */
2056 3469146 : vx = gvar(x);
2057 3469155 : vy = gvar(y);
2058 3469155 : if (vx != vy) /* x or y is treated as a scalar */
2059 2785619 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2060 2785619 : : mul_scal(y, x, ty);
2061 : /* vx = vy */
2062 2045731 : switch(tx)
2063 : {
2064 2045696 : case t_POL:
2065 : switch (ty)
2066 : {
2067 6741 : case t_SER:
2068 : {
2069 : long v;
2070 6741 : av = avma; v = RgX_valrem(x, &x);
2071 6741 : if (v == LONG_MAX) return gc_upto(av, Rg_get_0(x));
2072 6727 : v += valser(y); ly = lg(y);
2073 6727 : if (ly == 2) { set_avma(av); return zeroser(vx, v); }
2074 6496 : if (degpol(x))
2075 : {
2076 2030 : z = init_ser(ly, vx, v);
2077 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2078 2030 : return gc_GEN(av, fill_ser(z, x));
2079 : }
2080 : /* take advantage of x = c*t^v */
2081 4466 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2082 4466 : setvalser(y, v); return y;
2083 : }
2084 :
2085 2038955 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2086 : }
2087 0 : break;
2088 :
2089 35 : case t_SER:
2090 : switch (ty)
2091 : {
2092 35 : case t_RFRAC:
2093 35 : av = avma;
2094 35 : return gc_upto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2095 : }
2096 0 : break;
2097 : }
2098 0 : pari_err_TYPE2("*",x,y);
2099 : return NULL; /* LCOV_EXCL_LINE */
2100 : }
2101 :
2102 : /* return a nonnormalized result */
2103 : GEN
2104 112759 : sqr_ser_part(GEN x, long l1, long l2)
2105 : {
2106 : long i, j, l;
2107 : pari_sp av;
2108 : GEN Z, z, p1, p2;
2109 : long mi;
2110 112759 : if (l2 < l1) return zeroser(varn(x), 2*valser(x));
2111 112745 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2112 112745 : Z = cgetg(l2-l1+3, t_SER);
2113 112745 : Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
2114 112745 : z = Z + 2-l1;
2115 112745 : x += 2; mi = 0;
2116 425418 : for (i=0; i<l1; i++)
2117 : {
2118 312673 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2119 : }
2120 :
2121 720520 : for (i=l1; i<=l2; i++)
2122 : {
2123 607775 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2124 607775 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2125 1256988 : for (j=i-mi; j<=minss(l,mi); j++)
2126 649213 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2127 607775 : p1 = gshift(p1,1);
2128 607775 : if ((i&1) == 0 && p2[i>>1])
2129 93669 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2130 607775 : gel(z,i) = gc_upto(av,p1);
2131 : }
2132 112745 : return Z;
2133 : }
2134 :
2135 : /* (u + v X)^2 mod (X^2 + bX + c), b = 0 or -1
2136 : * = u^2 - c v^2 + (2uv - b v^2) X */
2137 : static GEN
2138 70 : sqrq(GEN x)
2139 : {
2140 70 : GEN T = gel(x,1), c = gel(T,2), b = gel(T,3);
2141 70 : GEN u = gel(x,2), v = gel(x,3), U, V, E, F, s, z;
2142 : pari_sp av, av2;
2143 :
2144 70 : z = cgetg(4, t_QUAD), gel(z,1) = ZX_copy(T); av = avma;
2145 70 : U = gsqr(u);
2146 70 : V = gsqr(v); s = gmul(c, V);
2147 70 : E = gmul(u, v);
2148 70 : F = signe(b)? gadd(E, V): E; /* u v - b v^2 */
2149 70 : av2 = avma;
2150 70 : gel(z,2) = gsub(U, s);
2151 70 : gel(z,3) = gadd(E, F);
2152 70 : gc_slice_unsafe(av, av2, z+2, 2); return z;
2153 : }
2154 :
2155 : GEN
2156 1183810047 : gsqr(GEN x)
2157 : {
2158 : long i, lx;
2159 : pari_sp av, tetpil;
2160 : GEN z, p1, p2, p3;
2161 :
2162 1183810047 : switch(typ(x))
2163 : {
2164 971113213 : case t_INT: return sqri(x);
2165 195633576 : case t_REAL: return sqrr(x);
2166 142408 : case t_INTMOD: { GEN X = gel(x,1);
2167 142408 : z = cgetg(3,t_INTMOD);
2168 142409 : gel(z,2) = gc_INT((pari_sp)z, remii(sqri(gel(x,2)), X));
2169 142409 : gel(z,1) = icopy(X); return z;
2170 : }
2171 4146814 : case t_FRAC: return sqrfrac(x);
2172 :
2173 8095929 : case t_COMPLEX:
2174 8095929 : if (isintzero(gel(x,1))) {
2175 240052 : av = avma;
2176 240052 : return gc_upto(av, gneg(gsqr(gel(x,2))));
2177 : }
2178 7855874 : z = cgetg(3,t_COMPLEX); av = avma;
2179 7855844 : p1 = gadd(gel(x,1),gel(x,2));
2180 7855616 : p2 = gsub(gel(x,1), gel(x,2));
2181 7855596 : p3 = gmul(gel(x,1),gel(x,2));
2182 7855787 : tetpil = avma;
2183 7855787 : gel(z,1) = gmul(p1,p2);
2184 7855812 : gel(z,2) = gshift(p3,1); gc_slice_unsafe(av,tetpil,z+1,2); return z;
2185 :
2186 4795 : case t_PADIC:
2187 : {
2188 4795 : GEN u = padic_u(x), p = padic_p(x), pd = padic_pd(x);
2189 4795 : long v2 = 2*valp(x), d = precp(x);
2190 4795 : if (!signe(u)) { x = gcopy(x); setvalp(x, v2); return x; }
2191 4739 : if (!absequaliu(p, 2))
2192 3206 : retmkpadic(Fp_sqr(u, pd), icopy(p), icopy(pd), v2, d);
2193 : /* p = 2*/
2194 1533 : if (d == 1) /* (1 + O(2))^2 = 1 + O(2^3) */
2195 7 : retmkpadic(gen_1, gen_2, utoipos(8), v2, 3);
2196 1526 : retmkpadic(remi2n(sqri(u), d + 1), gen_2, int2n(d + 1), v2, d + 1);
2197 : }
2198 70 : case t_QUAD: return sqrq(x);
2199 :
2200 419385 : case t_POLMOD: return sqr_polmod(gel(x,1), gel(x,2));
2201 :
2202 2967193 : case t_FFELT: return FF_sqr(x);
2203 :
2204 1278946 : case t_POL: return RgX_sqr(x);
2205 :
2206 35042 : case t_SER:
2207 35042 : lx = lg(x);
2208 35042 : if (ser_isexactzero(x)) {
2209 21 : GEN z = gcopy(x);
2210 21 : setvalser(z, 2*valser(x));
2211 21 : return z;
2212 : }
2213 35021 : if (lx < 40)
2214 34762 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2215 : else
2216 : {
2217 259 : pari_sp av = avma;
2218 259 : GEN z = init_ser(lx, varn(x), 2*valser(x));
2219 259 : x = ser2pol_i(x, lx);
2220 259 : x = RgXn_sqr(x, lx-2);
2221 259 : return gc_GEN(av, fill_ser(z,x));
2222 : }
2223 :
2224 63 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2225 63 : gel(z,1) = gsqr(gel(x,1));
2226 63 : gel(z,2) = gsqr(gel(x,2)); return z;
2227 :
2228 1169 : case t_MAT: return RgM_sqr(x);
2229 28 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
2230 : /* fall through handle extended t_QFB */
2231 36 : case t_QFB: return qfbsqr(x);
2232 658 : case t_VECSMALL:
2233 658 : z = cgetg_copy(x, &lx);
2234 16289 : for (i=1; i<lx; i++)
2235 : {
2236 15631 : long xi = x[i];
2237 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2238 15631 : z[i] = x[xi];
2239 : }
2240 658 : return z;
2241 : }
2242 7 : pari_err_TYPE2("*",x,x);
2243 : return NULL; /* LCOV_EXCL_LINE */
2244 : }
2245 :
2246 : /********************************************************************/
2247 : /** **/
2248 : /** DIVISION **/
2249 : /** **/
2250 : /********************************************************************/
2251 : static GEN
2252 32139 : div_rfrac_scal(GEN x, GEN y)
2253 : {
2254 32139 : pari_sp av = avma;
2255 32139 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2256 32139 : return gc_upto(av, gred_rfrac_simple(gel(x,1), d));
2257 : }
2258 : static GEN
2259 37767 : div_scal_rfrac(GEN x, GEN y)
2260 : {
2261 37767 : GEN y1 = gel(y,1), y2 = gel(y,2);
2262 37767 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2263 : {
2264 14 : if (degpol(y1))
2265 : {
2266 14 : pari_sp av = avma;
2267 14 : GEN _1 = Rg_get_1(x);
2268 14 : if (transtype(_1)) y1 = gmul(y1, _1);
2269 14 : return gc_upto(av, gred_rfrac_simple(gmul(x, y2), y1));
2270 : }
2271 0 : y1 = gel(y1,2);
2272 : }
2273 37753 : return RgX_Rg_mul(y2, gdiv(x,y1));
2274 : }
2275 : static GEN
2276 1187237 : div_rfrac(GEN x, GEN y)
2277 1187237 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2278 :
2279 : /* x != 0 */
2280 : static GEN
2281 1338593 : div_ser_scal(GEN x, GEN t)
2282 : {
2283 1338593 : if (ser_isexactzero(x))
2284 : {
2285 28 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gdiv(gel(x,2), t), varn(x), 1);
2286 28 : setvalser(z, valser(x)); return z;
2287 : }
2288 6202113 : pari_APPLY_ser(gdiv(gel(x,i), t));
2289 : }
2290 : GEN
2291 7658 : ser_normalize(GEN x)
2292 : {
2293 7658 : long i, lx = lg(x);
2294 : GEN c, z;
2295 7658 : if (lx == 2) return x;
2296 7658 : c = gel(x,2); if (gequal1(c)) return x;
2297 7581 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2298 108752 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2299 7581 : return z;
2300 : }
2301 :
2302 : /* y != 0 */
2303 : static GEN
2304 7124796 : div_T_scal(GEN x, GEN y, long tx) {
2305 7124796 : switch(tx)
2306 : {
2307 5757917 : case t_POL: return RgX_Rg_div(x, y);
2308 1338586 : case t_SER: return div_ser_scal(x, y);
2309 28296 : case t_RFRAC: return div_rfrac_scal(x,y);
2310 : }
2311 0 : pari_err_TYPE2("/",x,y);
2312 : return NULL; /* LCOV_EXCL_LINE */
2313 : }
2314 :
2315 : static GEN
2316 9257065 : div_scal_pol(GEN x, GEN y) {
2317 9257065 : long ly = lg(y);
2318 : pari_sp av;
2319 : GEN _1;
2320 9257065 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2321 9178167 : if (isrationalzero(x)) return zeropol(varn(y));
2322 7131383 : av = avma;
2323 7131383 : _1 = Rg_get_1(x); if (transtype(_1)) y = gmul(y, _1);
2324 7131383 : return gc_upto(av, gred_rfrac_simple(x,y));
2325 : }
2326 : static GEN
2327 18550 : div_scal_ser(GEN x, GEN y)
2328 : {
2329 18550 : pari_sp av = avma;
2330 18550 : GEN _1 = Rg_get_1(x);
2331 18550 : if (transtype(_1)) y = gmul(y, _1);
2332 18550 : return gc_upto(av, gmul(x, ser_inv(y)));
2333 : }
2334 : static GEN
2335 9264762 : div_scal_T(GEN x, GEN y, long ty) {
2336 9264762 : switch(ty)
2337 : {
2338 9209208 : case t_POL: return div_scal_pol(x, y);
2339 18543 : case t_SER: return div_scal_ser(x, y);
2340 37011 : case t_RFRAC: return div_scal_rfrac(x, y);
2341 : }
2342 0 : pari_err_TYPE2("/",x,y);
2343 : return NULL; /* LCOV_EXCL_LINE */
2344 : }
2345 :
2346 : /* assume tx = ty = t_SER, same variable vx */
2347 : static GEN
2348 771538 : div_ser(GEN x, GEN y, long vx)
2349 : {
2350 771538 : long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
2351 771538 : GEN y0 = y, z;
2352 771538 : pari_sp av = avma;
2353 :
2354 771538 : if (!signe(y)) pari_err_INV("div_ser", y);
2355 771531 : if (ser_isexactzero(x))
2356 : {
2357 59892 : if (lx == 2) return zeroser(vx, v);
2358 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2359 147 : setvalser(z, v); return z;
2360 : }
2361 711639 : if (lx < ly) ly = lx;
2362 711639 : y = ser2pol_i_normalize(y, ly, &e);
2363 711639 : if (e)
2364 : {
2365 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2366 7 : ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
2367 : }
2368 711639 : z = init_ser(ly, vx, v);
2369 711639 : if (ly == 3)
2370 : {
2371 12375 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2372 12375 : if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
2373 12375 : return gc_upto(av, z);
2374 : }
2375 699264 : x = ser2pol_i(x, ly);
2376 699264 : y = RgXn_div_i(x, y, ly-2);
2377 699264 : return gc_GEN(av, fill_ser(z,y));
2378 : }
2379 : /* x,y compatible PADIC */
2380 : static GEN
2381 13212327 : divpp(GEN x, GEN y)
2382 : {
2383 13212327 : GEN M, ux = padic_u(x), uy = padic_u(y), p = padic_p(x);
2384 13212327 : long a = precp(x), b = precp(y), v = valp(x) - valp(y);
2385 :
2386 13212327 : if (!signe(uy)) pari_err_INV("divpp",y);
2387 13212337 : if (!signe(ux)) return zeropadic(p, v);
2388 13212008 : if (a > b) M = padic_pd(y); else { M = padic_pd(x); b = a; }
2389 13212008 : retmkpadic(Fp_div(ux, uy, M), icopy(padic_p(x)), icopy(M), v, b);
2390 : }
2391 : static GEN
2392 36764 : div_polmod_same(GEN T, GEN x, GEN y)
2393 : {
2394 36764 : long v = varn(T);
2395 36764 : GEN a, z = cgetg(3, t_POLMOD);
2396 36764 : gel(z,1) = RgX_copy(T);
2397 36764 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2398 11102 : a = gdiv(x, y);
2399 25662 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2400 105 : {
2401 105 : pari_sp av = avma;
2402 105 : a = gc_upto(av, gmul(x, RgXQ_inv(y, T)));
2403 : }
2404 25557 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2405 5278 : {
2406 5278 : pari_sp av = avma;
2407 5278 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2408 5278 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2409 5278 : a = gc_upto(av, a);
2410 : }
2411 : else
2412 : {
2413 20279 : pari_sp av = avma;
2414 20279 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2415 20279 : a = gc_upto(av, a);
2416 : }
2417 36764 : gel(z,2) = a; return z;
2418 : }
2419 : GEN
2420 438953076 : gdiv(GEN x, GEN y)
2421 : {
2422 438953076 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2423 : pari_sp av, tetpil;
2424 : GEN z, p1;
2425 :
2426 438953076 : if (tx == ty) switch(tx)
2427 : {
2428 92086428 : case t_INT:
2429 92086428 : return Qdivii(x,y);
2430 :
2431 94229612 : case t_REAL: return divrr(x,y);
2432 25702 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2433 25702 : z = cgetg(3,t_INTMOD);
2434 25702 : if (X==Y || equalii(X,Y))
2435 25688 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2436 14 : gel(z,1) = gcdii(X,Y);
2437 14 : warn_coercion(X,Y,gel(z,1));
2438 14 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2439 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
2440 : }
2441 4043720 : case t_FRAC: {
2442 4043720 : GEN x1 = gel(x,1), x2 = gel(x,2);
2443 4043720 : GEN y1 = gel(y,1), y2 = gel(y,2);
2444 4043720 : z = cgetg(3, t_FRAC);
2445 4043720 : p1 = gcdii(x1, y1);
2446 4043717 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2447 4043716 : p1 = gcdii(x2, y2);
2448 4043712 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2449 4043711 : tetpil = avma;
2450 4043711 : gel(z,2) = mulii(x2,y1);
2451 4043715 : gel(z,1) = mulii(x1,y2);
2452 4043714 : normalize_frac(z);
2453 4043714 : fix_frac_if_int_GC(z,tetpil);
2454 4043718 : return z;
2455 : }
2456 9425515 : case t_COMPLEX:
2457 9425515 : if (isintzero(gel(y,1)))
2458 : {
2459 182777 : y = gel(y,2);
2460 182777 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2461 103943 : z = cgetg(3,t_COMPLEX);
2462 103943 : gel(z,1) = gdiv(gel(x,2), y);
2463 103943 : av = avma;
2464 103943 : gel(z,2) = gc_upto(av, gneg(gdiv(gel(x,1), y)));
2465 103943 : return z;
2466 : }
2467 9242738 : av = avma;
2468 9242738 : return gc_upto(av, gdiv(mulcc(x, conj_i(y)), cxnorm(y)));
2469 :
2470 587462 : case t_PADIC:
2471 587462 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("/",x,y);
2472 587462 : return divpp(x, y);
2473 :
2474 259 : case t_QUAD:
2475 259 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2476 245 : av = avma;
2477 245 : return gc_upto(av, gdiv(mulqq(x, conj_i(y)), quadnorm(y)));
2478 :
2479 246063 : case t_FFELT: return FF_div(x,y);
2480 :
2481 41139 : case t_POLMOD:
2482 41139 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2483 36764 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2484 : else {
2485 4375 : av = avma;
2486 4375 : z = gc_upto(av, gmul(x, ginv(y)));
2487 : }
2488 41139 : return z;
2489 :
2490 18571807 : case t_POL:
2491 18571807 : vx = varn(x);
2492 18571807 : vy = varn(y);
2493 18571807 : if (vx != vy) {
2494 154399 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2495 47857 : else return div_scal_pol(x, y);
2496 : }
2497 18417408 : if (!signe(y)) pari_err_INV("gdiv",y);
2498 18417408 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2499 18288403 : av = avma;
2500 18288403 : return gc_upto(av, gred_rfrac2(x,y));
2501 :
2502 771559 : case t_SER:
2503 771559 : vx = varn(x);
2504 771559 : vy = varn(y);
2505 771559 : if (vx != vy) {
2506 21 : if (varncmp(vx, vy) < 0)
2507 : {
2508 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2509 7 : return div_ser_scal(x, y);
2510 : }
2511 7 : return div_scal_ser(x, y);
2512 : }
2513 771538 : return div_ser(x, y, vx);
2514 1191696 : case t_RFRAC:
2515 1191696 : vx = varn(gel(x,2));
2516 1191696 : vy = varn(gel(y,2));
2517 1191696 : if (vx != vy) {
2518 4459 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2519 756 : else return div_scal_rfrac(x, y);
2520 : }
2521 1187237 : return div_rfrac(x,y);
2522 :
2523 21 : case t_VEC: /* handle extended t_QFB */
2524 21 : case t_QFB: av = avma; return gc_upto(av, qfbcomp(x, ginv(y)));
2525 :
2526 28 : case t_MAT:
2527 28 : p1 = RgM_div(x,y);
2528 28 : if (!p1) pari_err_INV("gdiv",y);
2529 21 : return p1;
2530 :
2531 0 : default: pari_err_TYPE2("/",x,y);
2532 : }
2533 :
2534 217736067 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2535 : {
2536 3804642 : long s = signe(x);
2537 3804642 : if (!s) {
2538 522269 : if (gequal0(y)) pari_err_INV("gdiv",y);
2539 522269 : switch (ty)
2540 : {
2541 519084 : default: return gen_0;
2542 70 : case t_INTMOD:
2543 70 : z = cgetg(3,t_INTMOD);
2544 70 : gel(z,1) = icopy(gel(y,1));
2545 70 : gel(z,2) = gen_0; return z;
2546 3115 : case t_FFELT: return FF_zero(y);
2547 : }
2548 : }
2549 3282373 : if (is_pm1(x)) {
2550 1323615 : if (s > 0) return ginv(y);
2551 228729 : av = avma; return gc_upto(av, ginv(gneg(y)));
2552 : }
2553 1958757 : switch(ty)
2554 : {
2555 644380 : case t_REAL: return divir(x,y);
2556 21 : case t_INTMOD:
2557 21 : z = cgetg(3, t_INTMOD);
2558 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2559 792515 : case t_FRAC:
2560 792515 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2561 792515 : if (equali1(p1))
2562 : {
2563 271589 : set_avma((pari_sp)z);
2564 271589 : gel(z,2) = icopy(gel(y,1));
2565 271589 : gel(z,1) = mulii(gel(y,2), x);
2566 271589 : normalize_frac(z);
2567 271589 : fix_frac_if_int(z);
2568 : }
2569 : else
2570 : {
2571 520926 : x = diviiexact(x,p1); tetpil = avma;
2572 520926 : gel(z,2) = diviiexact(gel(y,1), p1);
2573 520926 : gel(z,1) = mulii(gel(y,2), x);
2574 520926 : normalize_frac(z);
2575 520926 : fix_frac_if_int_GC(z,tetpil);
2576 : }
2577 792515 : return z;
2578 :
2579 420 : case t_FFELT: return Z_FF_div(x,y);
2580 519694 : case t_COMPLEX: return divRc(x,y);
2581 1505 : case t_PADIC: return divTp(x, y);
2582 231 : case t_QUAD: return divRq(x,y);
2583 : }
2584 : }
2585 213931415 : if (gequal0(y))
2586 : {
2587 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2588 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2589 : }
2590 :
2591 213939615 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2592 : {
2593 173327683 : case t_REAL:
2594 173327683 : switch(ty)
2595 : {
2596 171046303 : case t_INT: return divri(x,y);
2597 519535 : case t_FRAC:
2598 519535 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2599 519535 : return gc_leaf(av, z);
2600 1761891 : case t_COMPLEX: return divRc(x, y);
2601 42 : case t_QUAD: return divfq(x, y, realprec(x));
2602 18 : default: pari_err_TYPE2("/",x,y);
2603 : }
2604 :
2605 273 : case t_INTMOD:
2606 : switch(ty)
2607 : {
2608 189 : case t_INT:
2609 189 : z = cgetg(3, t_INTMOD);
2610 189 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2611 28 : case t_FRAC: { GEN X = gel(x,1);
2612 28 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2613 28 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2614 : }
2615 7 : case t_FFELT:
2616 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2617 0 : pari_err_OP("/",x,y);
2618 7 : return Z_FF_div(gel(x,2),y);
2619 :
2620 0 : case t_COMPLEX: return divRc(x,y);
2621 14 : case t_QUAD: return divRq(x,y);
2622 :
2623 7 : case t_PADIC: { GEN X = gel(x,1);
2624 7 : z = cgetg(3, t_INTMOD);
2625 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2626 : }
2627 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2628 : }
2629 :
2630 : case t_FRAC:
2631 : switch(ty)
2632 : {
2633 2102453 : case t_INT: z = cgetg(3, t_FRAC);
2634 2102453 : p1 = gcdii(y,gel(x,1));
2635 2102451 : if (equali1(p1))
2636 : {
2637 832588 : set_avma((pari_sp)z); tetpil = 0;
2638 832588 : gel(z,1) = icopy(gel(x,1));
2639 : }
2640 : else
2641 : {
2642 1269865 : y = diviiexact(y,p1); tetpil = avma;
2643 1269865 : gel(z,1) = diviiexact(gel(x,1), p1);
2644 : }
2645 2102452 : gel(z,2) = mulii(gel(x,2),y);
2646 2102451 : normalize_frac(z);
2647 2102451 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2648 2102453 : return z;
2649 :
2650 59572 : case t_REAL:
2651 59572 : av = avma;
2652 59572 : return gc_leaf(av, divir(gel(x,1), mulri(y,gel(x,2))));
2653 :
2654 7 : case t_INTMOD: { GEN Y = gel(y,1);
2655 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2656 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2657 : }
2658 :
2659 28 : case t_FFELT: av=avma;
2660 28 : return gc_upto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2661 :
2662 868 : case t_COMPLEX: return divRc(x, y);
2663 :
2664 2141 : case t_PADIC:
2665 2141 : if (!signe(gel(x,1))) return gen_0;
2666 2141 : return divTp(x, y);
2667 :
2668 42 : case t_QUAD: return divRq(x, y);
2669 : }
2670 :
2671 : case t_FFELT:
2672 161 : switch (ty)
2673 : {
2674 77 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2675 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2676 7 : case t_INTMOD:
2677 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2678 0 : pari_err_OP("/",x,y);
2679 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2680 49 : default:
2681 49 : pari_err_TYPE2("/",x,y);
2682 : }
2683 0 : break;
2684 :
2685 13735547 : case t_COMPLEX:
2686 : switch(ty)
2687 : {
2688 13735546 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2689 0 : case t_INTMOD: return mulRc(ginv(y), x);
2690 0 : case t_PADIC:
2691 0 : return Zp_nosquare_m1(padic_p(y))? divcR(x,y): divTp(x, y);
2692 0 : case t_QUAD:
2693 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2694 0 : return divfq(x, y, lx);
2695 : }
2696 :
2697 : case t_PADIC:
2698 : switch(ty)
2699 : {
2700 1205008 : case t_INT: case t_FRAC: { GEN p = padic_p(x);
2701 1204917 : return signe(padic_u(x))? divpT(x, y)
2702 2409925 : : zeropadic(p, valp(x) - Q_pval(y,p));
2703 : }
2704 7 : case t_INTMOD: { GEN Y = gel(y,1);
2705 7 : z = cgetg(3, t_INTMOD);
2706 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2707 : }
2708 0 : case t_COMPLEX: return divRc(x,y);
2709 14 : case t_QUAD: return divRq(x,y);
2710 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2711 : }
2712 :
2713 : case t_QUAD:
2714 : switch (ty)
2715 : {
2716 1204 : case t_INT: case t_INTMOD: case t_FRAC:
2717 1204 : z = cgetg(4,t_QUAD);
2718 1204 : gel(z,1) = ZX_copy(gel(x,1));
2719 1204 : gel(z,2) = gdiv(gel(x,2), y);
2720 1204 : gel(z,3) = gdiv(gel(x,3), y); return z;
2721 63 : case t_REAL: return divqf(x, y, realprec(y));
2722 28 : case t_PADIC: return divTp(x, y);
2723 0 : case t_COMPLEX:
2724 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2725 0 : return divqf(x, y, ly);
2726 : }
2727 : }
2728 23504398 : switch(ty) {
2729 582625 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2730 582625 : av = avma; return gc_upto(av, gmul(x, ginv(y)));
2731 42 : case t_MAT:
2732 42 : av = avma; p1 = RgM_inv(y);
2733 35 : if (!p1) pari_err_INV("gdiv",y);
2734 28 : return gc_upto(av, gmul(x, p1));
2735 0 : case t_VEC: case t_COL:
2736 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2737 0 : pari_err_TYPE2("/",x,y);
2738 : }
2739 22921731 : switch(tx) {
2740 4422483 : case t_VEC: case t_COL: case t_MAT:
2741 14369509 : pari_APPLY_same(gdiv(gel(x,i),y));
2742 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2743 0 : pari_err_TYPE2("/",x,y);
2744 : }
2745 :
2746 18499248 : vy = gvar(y);
2747 18500148 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2748 218130 : vx = varn(X);
2749 218130 : if (vx != vy) {
2750 217360 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2751 216947 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2752 : }
2753 : /* y is POL, SER or RFRAC */
2754 770 : av = avma;
2755 770 : switch(ty)
2756 : {
2757 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2758 770 : default: y = ginvmod(gmod(y,X), X);
2759 : }
2760 763 : return gc_upto(av, mul_polmod_same(X, gel(x,2), y));
2761 : }
2762 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2763 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2764 : * variable NO_VARIABLE, then the operation is incorrect */
2765 18282018 : vx = gvar(x);
2766 18282008 : if (vx != vy) { /* includes cases where one is scalar */
2767 16389149 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2768 9264348 : else return div_scal_T(x, y, ty);
2769 : }
2770 1892859 : switch(tx)
2771 : {
2772 1046171 : case t_POL:
2773 : switch(ty)
2774 : {
2775 28 : case t_SER:
2776 : {
2777 28 : GEN y0 = y;
2778 : long v;
2779 28 : av = avma; v = RgX_valrem(x, &x);
2780 28 : if (v == LONG_MAX) return gc_upto(av, Rg_get_0(x));
2781 14 : v -= valser(y); ly = lg(y); /* > 2 */
2782 14 : y = ser2pol_i_normalize(y, ly, &i);
2783 14 : if (i)
2784 : {
2785 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2786 7 : ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
2787 : }
2788 7 : z = init_ser(ly, vx, v);
2789 7 : return gc_GEN(av, fill_ser(z, RgXn_div(x, y, ly-2)));
2790 : }
2791 :
2792 1046143 : case t_RFRAC:
2793 : {
2794 1046143 : GEN y1 = gel(y,1), y2 = gel(y,2);
2795 1046143 : if (typ(y1) == t_POL && varn(y1) == vx)
2796 1143 : return mul_rfrac_scal(y2, y1, x);
2797 1045000 : av = avma;
2798 1045000 : return gc_upto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2799 : }
2800 : }
2801 0 : break;
2802 :
2803 584324 : case t_SER:
2804 : switch(ty)
2805 : {
2806 584310 : case t_POL:
2807 : {
2808 584310 : long v = valser(x);
2809 584310 : lx = lg(x);
2810 584310 : if (lx == 2) return zeroser(vx, v - RgX_val(y));
2811 584205 : av = avma;
2812 584205 : x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
2813 584205 : z = init_ser(lx, vx, v);
2814 584205 : if (!signe(x)) setsigne(z,0);
2815 584205 : return gc_GEN(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
2816 : }
2817 14 : case t_RFRAC:
2818 14 : av = avma;
2819 14 : return gc_upto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2820 : }
2821 0 : break;
2822 :
2823 262344 : case t_RFRAC:
2824 : switch(ty)
2825 : {
2826 262344 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2827 0 : case t_SER:
2828 0 : av = avma;
2829 0 : return gc_upto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2830 : }
2831 0 : break;
2832 : }
2833 20 : pari_err_TYPE2("/",x,y);
2834 : return NULL; /* LCOV_EXCL_LINE */
2835 : }
2836 :
2837 : /********************************************************************/
2838 : /** **/
2839 : /** SIMPLE MULTIPLICATION **/
2840 : /** **/
2841 : /********************************************************************/
2842 : GEN
2843 236592560 : gmulsg(long s, GEN x)
2844 : {
2845 : pari_sp av;
2846 : long i;
2847 : GEN z;
2848 :
2849 236592560 : switch(typ(x))
2850 : {
2851 142462885 : case t_INT: return mulsi(s,x);
2852 70384092 : case t_REAL: return s? mulsr(s,x): gen_0; /* gmul semantic */
2853 174308 : case t_INTMOD: { GEN p = gel(x,1);
2854 174308 : z = cgetg(3,t_INTMOD);
2855 174308 : gel(z,2) = gc_INT((pari_sp)z, modii(mulsi(s,gel(x,2)), p));
2856 174307 : gel(z,1) = icopy(p); return z;
2857 : }
2858 548125 : case t_FFELT: return FF_Z_mul(x,stoi(s));
2859 6159113 : case t_FRAC:
2860 6159113 : if (!s) return gen_0;
2861 6148417 : z = cgetg(3,t_FRAC);
2862 6148417 : i = labs(s); i = ugcdiu(gel(x,2), i);
2863 6148417 : if (i == 1)
2864 : {
2865 2259746 : gel(z,2) = icopy(gel(x,2));
2866 2259746 : gel(z,1) = mulis(gel(x,1), s);
2867 : }
2868 : else
2869 : {
2870 3888671 : gel(z,2) = diviuexact(gel(x,2), (ulong)i);
2871 3888671 : gel(z,1) = mulis(gel(x,1), s/i);
2872 3888671 : fix_frac_if_int(z);
2873 : }
2874 6148417 : return z;
2875 :
2876 14537397 : case t_COMPLEX:
2877 14537397 : if (!s) return gen_0;
2878 13411593 : z = cgetg(3, t_COMPLEX);
2879 13411546 : gel(z,1) = gmulsg(s,gel(x,1));
2880 13409973 : gel(z,2) = gmulsg(s,gel(x,2)); return z;
2881 :
2882 1449 : case t_PADIC:
2883 1449 : if (!s) return gen_0;
2884 1449 : av = avma; return gc_upto(av, mulpp(cvtop2(stoi(s),x), x));
2885 :
2886 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2887 7 : gel(z,1) = ZX_copy(gel(x,1));
2888 7 : gel(z,2) = gmulsg(s,gel(x,2));
2889 7 : gel(z,3) = gmulsg(s,gel(x,3)); return z;
2890 :
2891 208516 : case t_POLMOD:
2892 208516 : retmkpolmod(gmulsg(s,gel(x,2)), RgX_copy(gel(x,1)));
2893 :
2894 821879 : case t_POL:
2895 821879 : if (!signe(x)) return RgX_copy(x);
2896 801495 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2897 3329122 : pari_APPLY_pol(gmulsg(s,gel(x,i)));
2898 :
2899 182 : case t_SER:
2900 182 : if (ser_isexactzero(x)) return gcopy(x);
2901 182 : if (!s) return Rg_get_0(x);
2902 3864 : pari_APPLY_ser(gmulsg(s,gel(x,i)));
2903 :
2904 0 : case t_RFRAC:
2905 0 : if (!s) return zeropol(varn(gel(x,2)));
2906 0 : if (s == 1) return gcopy(x);
2907 0 : if (s == -1) return gneg(x);
2908 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), stoi(s));
2909 :
2910 1304585 : case t_VEC: case t_COL: case t_MAT:
2911 4109873 : pari_APPLY_same(gmulsg(s,gel(x,i)));
2912 : }
2913 0 : pari_err_TYPE("gmulsg",x);
2914 : return NULL; /* LCOV_EXCL_LINE */
2915 : }
2916 :
2917 : GEN
2918 282636557 : gmulug(ulong s, GEN x)
2919 : {
2920 : pari_sp av;
2921 : long i;
2922 : GEN z;
2923 :
2924 282636557 : switch(typ(x))
2925 : {
2926 280321653 : case t_INT: return mului(s,x);
2927 1064579 : case t_REAL: return s? mulur(s,x): gen_0; /* gmul semantic */
2928 364 : case t_INTMOD: { GEN p = gel(x,1);
2929 364 : z = cgetg(3,t_INTMOD);
2930 364 : gel(z,2) = gc_INT((pari_sp)z, modii(mului(s,gel(x,2)), p));
2931 364 : gel(z,1) = icopy(p); return z;
2932 : }
2933 413 : case t_FFELT: return FF_Z_mul(x,utoi(s));
2934 981585 : case t_FRAC:
2935 981585 : if (!s) return gen_0;
2936 981571 : z = cgetg(3,t_FRAC);
2937 981571 : i = ugcdiu(gel(x,2), s);
2938 981571 : if (i == 1)
2939 : {
2940 822559 : gel(z,2) = icopy(gel(x,2));
2941 822559 : gel(z,1) = muliu(gel(x,1), s);
2942 : }
2943 : else
2944 : {
2945 159012 : gel(z,2) = diviuexact(gel(x,2), i);
2946 159012 : gel(z,1) = muliu(gel(x,1), s/i);
2947 159012 : fix_frac_if_int(z);
2948 : }
2949 981571 : return z;
2950 :
2951 35945 : case t_COMPLEX:
2952 35945 : if (!s) return gen_0;
2953 35945 : z = cgetg(3, t_COMPLEX);
2954 35945 : gel(z,1) = gmulug(s,gel(x,1));
2955 35945 : gel(z,2) = gmulug(s,gel(x,2)); return z;
2956 :
2957 7342 : case t_PADIC:
2958 7342 : if (!s) return gen_0;
2959 7342 : av = avma; return gc_upto(av, mulpp(cvtop2(utoi(s),x), x));
2960 :
2961 0 : case t_QUAD: z = cgetg(4, t_QUAD);
2962 0 : gel(z,1) = ZX_copy(gel(x,1));
2963 0 : gel(z,2) = gmulug(s,gel(x,2));
2964 0 : gel(z,3) = gmulug(s,gel(x,3)); return z;
2965 :
2966 6363 : case t_POLMOD:
2967 6363 : retmkpolmod(gmulug(s,gel(x,2)), RgX_copy(gel(x,1)));
2968 :
2969 200697 : case t_POL:
2970 200697 : if (!signe(x)) return RgX_copy(x);
2971 199920 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2972 621229 : pari_APPLY_pol(gmulug(s,gel(x,i)));
2973 :
2974 0 : case t_SER:
2975 0 : if (ser_isexactzero(x)) return gcopy(x);
2976 0 : if (!s) return Rg_get_0(x);
2977 0 : pari_APPLY_ser(gmulug(s,gel(x,i)));
2978 :
2979 0 : case t_RFRAC:
2980 0 : if (!s) return zeropol(varn(gel(x,2)));
2981 0 : if (s == 1) return gcopy(x);
2982 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), utoi(s));
2983 :
2984 17618 : case t_VEC: case t_COL: case t_MAT:
2985 96479 : pari_APPLY_same(gmulug(s,gel(x,i)));
2986 : }
2987 0 : pari_err_TYPE("gmulsg",x);
2988 : return NULL; /* LCOV_EXCL_LINE */
2989 : }
2990 :
2991 : /********************************************************************/
2992 : /** **/
2993 : /** SIMPLE DIVISION **/
2994 : /** **/
2995 : /********************************************************************/
2996 :
2997 : GEN
2998 12995248 : gdivgs(GEN x, long s)
2999 : {
3000 12995248 : long tx = typ(x), i;
3001 : pari_sp av;
3002 : GEN z;
3003 :
3004 12995248 : if (!s)
3005 : {
3006 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3007 0 : pari_err_INV("gdivgs",gen_0);
3008 : }
3009 12995261 : switch(tx)
3010 : {
3011 1592020 : case t_INT: return Qdivis(x, s);
3012 8553360 : case t_REAL: return divrs(x,s);
3013 :
3014 357 : case t_INTMOD:
3015 357 : z = cgetg(3, t_INTMOD);
3016 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3017 :
3018 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3019 :
3020 556873 : case t_FRAC: z = cgetg(3, t_FRAC);
3021 556873 : i = ugcdiu(gel(x,1), labs(s));
3022 556873 : if (i == 1)
3023 : {
3024 413436 : gel(z,2) = mulsi(s, gel(x,2));
3025 413436 : gel(z,1) = icopy(gel(x,1));
3026 : }
3027 : else
3028 : {
3029 143437 : gel(z,2) = mulsi(s/i, gel(x,2));
3030 143437 : gel(z,1) = divis(gel(x,1), i);
3031 : }
3032 556873 : normalize_frac(z);
3033 556873 : fix_frac_if_int(z); return z;
3034 :
3035 1809612 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3036 1809612 : gel(z,1) = gdivgs(gel(x,1),s);
3037 1809608 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3038 :
3039 133 : case t_PADIC: /* divpT */
3040 : {
3041 133 : GEN p = padic_p(x);
3042 133 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3043 133 : av = avma;
3044 133 : return gc_upto(av, divpp(x, cvtop2(stoi(s),x)));
3045 : }
3046 :
3047 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3048 28 : gel(z,1) = ZX_copy(gel(x,1));
3049 28 : gel(z,2) = gdivgs(gel(x,2),s);
3050 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3051 :
3052 37282 : case t_POLMOD:
3053 37282 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3054 :
3055 91 : case t_RFRAC:
3056 91 : if (s == 1) return gcopy(x);
3057 84 : else if (s == -1) return gneg(x);
3058 84 : return div_rfrac_scal(x, stoi(s));
3059 :
3060 156754 : case t_POL: pari_APPLY_pol_normalized(gdivgs(gel(x,i),s));
3061 0 : case t_SER: pari_APPLY_ser_normalized(gdivgs(gel(x,i),s));
3062 407593 : case t_VEC:
3063 : case t_COL:
3064 894977 : case t_MAT: pari_APPLY_same(gdivgs(gel(x,i),s));
3065 : }
3066 0 : pari_err_TYPE2("/",x, stoi(s));
3067 : return NULL; /* LCOV_EXCL_LINE */
3068 : }
3069 :
3070 : GEN
3071 60356051 : gdivgu(GEN x, ulong s)
3072 : {
3073 60356051 : long tx = typ(x), i;
3074 : pari_sp av;
3075 : GEN z;
3076 :
3077 60356051 : if (!s)
3078 : {
3079 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3080 0 : pari_err_INV("gdivgu",gen_0);
3081 : }
3082 60356003 : switch(tx)
3083 : {
3084 23521084 : case t_INT: return Qdiviu(x, s);
3085 12899467 : case t_REAL: return divru(x,s);
3086 :
3087 210315 : case t_INTMOD:
3088 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3089 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3090 :
3091 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3092 :
3093 727752 : case t_FRAC: z = cgetg(3, t_FRAC);
3094 727752 : i = ugcdiu(gel(x,1), s);
3095 727752 : if (i == 1)
3096 : {
3097 534167 : gel(z,2) = mului(s, gel(x,2));
3098 534167 : gel(z,1) = icopy(gel(x,1));
3099 : }
3100 : else
3101 : {
3102 193585 : gel(z,2) = mului(s/i, gel(x,2));
3103 193585 : gel(z,1) = divis(gel(x,1), i);
3104 : }
3105 727752 : normalize_frac(z);
3106 727752 : fix_frac_if_int(z); return z;
3107 :
3108 11347210 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3109 11347210 : gel(z,1) = gdivgu(gel(x,1),s);
3110 11347210 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3111 :
3112 11551237 : case t_PADIC: /* divpT */
3113 : {
3114 11551237 : GEN p = padic_p(x);
3115 11551237 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3116 11415703 : av = avma;
3117 11415703 : return gc_upto(av, divpp(x, cvtop2(utoi(s),x)));
3118 : }
3119 :
3120 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3121 0 : gel(z,1) = ZX_copy(gel(x,1));
3122 0 : gel(z,2) = gdivgu(gel(x,2),s);
3123 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3124 :
3125 1456 : case t_POLMOD:
3126 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3127 :
3128 56 : case t_RFRAC:
3129 56 : if (s == 1) return gcopy(x);
3130 56 : return div_rfrac_scal(x, utoi(s));
3131 :
3132 450618 : case t_POL: pari_APPLY_pol_normalized(gdivgu(gel(x,i),s));
3133 16786 : case t_SER: pari_APPLY_ser_normalized(gdivgu(gel(x,i),s));
3134 336 : case t_VEC:
3135 : case t_COL:
3136 1162 : case t_MAT: pari_APPLY_same(gdivgu(gel(x,i),s));
3137 : }
3138 0 : pari_err_TYPE2("/",x, utoi(s));
3139 : return NULL; /* LCOV_EXCL_LINE */
3140 : }
3141 :
3142 : /* x / (i*(i+1)) */
3143 : GEN
3144 224428476 : divrunextu(GEN x, ulong i)
3145 : {
3146 224428476 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3147 0 : return divri(x, muluu(i , i+1));
3148 : else
3149 224428476 : return divru(x, i*(i+1));
3150 : }
3151 : /* x / (i*(i+1)) */
3152 : GEN
3153 1583873 : gdivgunextu(GEN x, ulong i)
3154 : {
3155 1583873 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3156 0 : return gdivgu(x, i*(i+1));
3157 : else
3158 1583873 : return gdiv(x, muluu(i, i+1));
3159 : }
3160 :
3161 : /* True shift (exact multiplication by 2^n) */
3162 : GEN
3163 143493455 : gmul2n(GEN x, long n)
3164 : {
3165 : GEN z, a, b;
3166 : long k, l;
3167 :
3168 143493455 : switch(typ(x))
3169 : {
3170 43038119 : case t_INT:
3171 43038119 : if (n>=0) return shifti(x,n);
3172 7634457 : if (!signe(x)) return gen_0;
3173 5080124 : l = vali(x); n = -n;
3174 5080193 : if (n<=l) return shifti(x,-n);
3175 348553 : z = cgetg(3,t_FRAC);
3176 348553 : gel(z,1) = shifti(x,-l);
3177 348553 : gel(z,2) = int2n(n-l); return z;
3178 :
3179 64249337 : case t_REAL:
3180 64249337 : return shiftr(x,n);
3181 :
3182 180937 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3183 180937 : z = cgetg(3,t_INTMOD);
3184 180937 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3185 82776 : gel(z,2) = gc_INT((pari_sp)z, modii(shifti(a,n), b));
3186 82773 : gel(z,1) = icopy(b); return z;
3187 :
3188 217433 : case t_FFELT: return FF_mul2n(x,n);
3189 :
3190 1128264 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3191 1128264 : l = vali(a);
3192 1128264 : k = vali(b);
3193 1128264 : if (n+l >= k)
3194 : {
3195 322881 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3196 229378 : l = n-k; k = -k;
3197 : }
3198 : else
3199 : {
3200 805383 : k = -(l+n); l = -l;
3201 : }
3202 1034761 : z = cgetg(3,t_FRAC);
3203 1034761 : gel(z,1) = shifti(a,l);
3204 1034761 : gel(z,2) = shifti(b,k); return z;
3205 :
3206 11059786 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3207 11059788 : gel(z,1) = gmul2n(gel(x,1),n);
3208 11059787 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3209 :
3210 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3211 105 : gel(z,1) = ZX_copy(gel(x,1));
3212 105 : gel(z,2) = gmul2n(gel(x,2),n);
3213 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3214 :
3215 177009 : case t_POLMOD:
3216 177009 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3217 :
3218 6752887 : case t_POL:
3219 24035352 : pari_APPLY_pol(gmul2n(gel(x,i),n));
3220 102875 : case t_SER:
3221 102875 : if (ser_isexactzero(x)) return gcopy(x);
3222 634312 : pari_APPLY_ser(gmul2n(gel(x,i),n));
3223 16580929 : case t_VEC: case t_COL: case t_MAT:
3224 61002045 : pari_APPLY_same(gmul2n(gel(x,i),n));
3225 :
3226 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3227 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3228 :
3229 7644 : case t_PADIC: /* int2n wrong if n < 0 */
3230 7644 : return gmul(gmul2n(gen_1,n),x);
3231 : }
3232 0 : pari_err_TYPE("gmul2n",x);
3233 : return NULL; /* LCOV_EXCL_LINE */
3234 : }
3235 :
3236 : /*******************************************************************/
3237 : /* */
3238 : /* INVERSE */
3239 : /* */
3240 : /*******************************************************************/
3241 : static GEN
3242 189627 : inv_polmod(GEN T, GEN x)
3243 : {
3244 189627 : GEN z = cgetg(3,t_POLMOD), a;
3245 189627 : gel(z,1) = RgX_copy(T);
3246 189627 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3247 64917 : a = ginv(x);
3248 : else
3249 : {
3250 124710 : if (lg(T) == 5) /* quadratic fields */
3251 13125 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3252 : else
3253 111585 : a = RgXQ_inv(x, T);
3254 : }
3255 189627 : gel(z,2) = a; return z;
3256 : }
3257 :
3258 : GEN
3259 36580876 : ginv(GEN x)
3260 : {
3261 : long s;
3262 : pari_sp av;
3263 : GEN z, y;
3264 :
3265 36580876 : switch(typ(x))
3266 : {
3267 9506422 : case t_INT:
3268 9506422 : if (is_pm1(x)) return icopy(x);
3269 7958878 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3270 7958864 : z = cgetg(3,t_FRAC);
3271 7965383 : gel(z,1) = s<0? gen_m1: gen_1;
3272 7965383 : gel(z,2) = absi(x); return z;
3273 :
3274 10867884 : case t_REAL: return invr(x);
3275 :
3276 5733 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3277 5733 : gel(z,1) = icopy(gel(x,1));
3278 5733 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3279 :
3280 533036 : case t_FRAC: {
3281 533036 : GEN a = gel(x,1), b = gel(x,2);
3282 533036 : s = signe(a);
3283 533036 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3284 247495 : z = cgetg(3,t_FRAC);
3285 247495 : gel(z,1) = icopy(b);
3286 247495 : gel(z,2) = icopy(a);
3287 247495 : normalize_frac(z); return z;
3288 : }
3289 9617346 : case t_COMPLEX:
3290 9617346 : av = avma;
3291 9617346 : return gc_upto(av, divcR(conj_i(x), cxnorm(x)));
3292 :
3293 273 : case t_QUAD:
3294 273 : av = avma;
3295 273 : return gc_upto(av, gdiv(conj_i(x), quadnorm(x)));
3296 :
3297 358358 : case t_PADIC:
3298 : {
3299 358358 : GEN p = padic_p(x), pd = padic_pd(x), u = padic_u(x);
3300 358358 : long d = precp(x);
3301 358358 : if (!signe(u)) pari_err_INV("ginv",x);
3302 358351 : retmkpadic(Zp_inv(u, p, d), icopy(p), icopy(pd), -valp(x), d);
3303 : }
3304 :
3305 189627 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3306 14069 : case t_FFELT: return FF_inv(x);
3307 5322037 : case t_POL: return gred_rfrac_simple(gen_1,x);
3308 26369 : case t_SER: return ser_inv(x);
3309 2933 : case t_RFRAC:
3310 : {
3311 2933 : GEN n = gel(x,1), d = gel(x,2);
3312 2933 : pari_sp av = avma, ltop;
3313 2933 : if (gequal0(n)) pari_err_INV("ginv",x);
3314 :
3315 2933 : n = simplify_shallow(n);
3316 2933 : if (typ(n) != t_POL || varn(n) != varn(d))
3317 : {
3318 2933 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3319 679 : ltop = avma;
3320 679 : z = RgX_Rg_div(d,n);
3321 : } else {
3322 0 : ltop = avma;
3323 0 : z = cgetg(3,t_RFRAC);
3324 0 : gel(z,1) = RgX_copy(d);
3325 0 : gel(z,2) = RgX_copy(n);
3326 : }
3327 679 : stackdummy(av, ltop);
3328 679 : return z;
3329 : }
3330 :
3331 21 : case t_VEC: if (!is_ext_qfr(x)) break;
3332 : case t_QFB:
3333 98910 : return qfbpow(x, gen_m1);
3334 38932 : case t_MAT:
3335 38932 : y = RgM_inv(x);
3336 38925 : if (!y) pari_err_INV("ginv",x);
3337 38855 : return y;
3338 28 : case t_VECSMALL:
3339 : {
3340 28 : long i, lx = lg(x)-1;
3341 28 : y = zero_zv(lx);
3342 112 : for (i=1; i<=lx; i++)
3343 : {
3344 84 : long xi = x[i];
3345 84 : if (xi<1 || xi>lx || y[xi])
3346 0 : pari_err_TYPE("ginv [not a permutation]", x);
3347 84 : y[xi] = i;
3348 : }
3349 28 : return y;
3350 : }
3351 : }
3352 6 : pari_err_TYPE("inverse",x);
3353 : return NULL; /* LCOV_EXCL_LINE */
3354 : }
|
