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 306668483 : addRc(GEN x, GEN y) {
61 306668483 : GEN z = cgetg(3,t_COMPLEX);
62 306648396 : gel(z,1) = gadd(x,gel(y,1));
63 306629398 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 353776410 : mulRc(GEN x, GEN y)
67 : {
68 353776410 : GEN a = gel(y,1), b = gel(y,2), z = cgetg(3,t_COMPLEX);
69 353780971 : gel(z,1) = (typ(x) != t_INTMOD && isintzero(a))? gen_0: gmul(x, a);
70 354251422 : gel(z,2) = gmul(x, b); return z;
71 : }
72 : static GEN
73 2365558 : divRc(GEN x, GEN y)
74 : {
75 2365558 : GEN a = gel(y,1), b = gel(y,2), z = cgetg(3,t_COMPLEX);
76 2365557 : pari_sp av = avma, av2;
77 2365557 : if (isintzero(a) && typ(x) != t_INTMOD)
78 : {
79 79282 : gel(z,1) = gen_0;
80 79282 : gel(z,2) = gc_upto(av, gdiv(x, gneg(b)));
81 : }
82 : else
83 : {
84 2286275 : GEN t = gdiv(x, cxnorm(y)), mb = gneg(b);
85 2286272 : av2 = avma;
86 2286272 : gel(z,1) = gmul(t, a);
87 2286262 : gel(z,2) = gmul(t, mb);
88 2286266 : gc_slice_unsafe(av, av2, z+1, 2);
89 : }
90 2365562 : 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 23399075 : divcR(GEN x, GEN y) {
100 23399075 : GEN z = cgetg(3,t_COMPLEX);
101 23399003 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
102 23398582 : 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 37790594 : mulrfrac(GEN x, GEN y)
126 : {
127 : pari_sp av;
128 37790594 : GEN z, a = gel(y,1), b = gel(y,2);
129 37790594 : if (is_pm1(a)) /* frequent special case */
130 : {
131 12769990 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
132 12769586 : return z;
133 : }
134 25020386 : av = avma;
135 25020386 : 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 4984905 : mulTp(GEN x, GEN y) { pari_sp av = avma;
163 4984905 : 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 3344850 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 3344850 : if (lgefint(X) == 3) {
182 3021727 : ulong u = Fl_add(itou(x),itou(y), X[2]);
183 3021727 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
184 : }
185 : else {
186 323123 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
187 323123 : gel(z,2) = gc_INT((pari_sp)z, u);
188 : }
189 3344850 : 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 2629232 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
206 2629232 : if (lgefint(X) == 3) {
207 2487916 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
208 2487916 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 141316 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x,y), X) );
212 2629229 : gel(z,1) = icopy(X); return z;
213 : }
214 : /* cf add_intmod_same */
215 : static GEN
216 341930 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
217 : {
218 341930 : if (lgefint(X) == 3) {
219 319390 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
220 319327 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
221 : }
222 : else
223 22540 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
224 341867 : 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 9959372 : rfrac_denom_mul_scal(GEN d, GEN y)
237 : {
238 9959372 : GEN D = RgX_Rg_mul(d, y);
239 9959372 : 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 9959372 : return D;
245 : }
246 :
247 : static int
248 57792533 : Leading_is_neg(GEN x)
249 : {
250 122218885 : while (typ(x) == t_POL) x = leading_coeff(x);
251 57792533 : return is_real_t(typ(x))? gsigne(x) < 0: 0;
252 : }
253 :
254 : static int
255 154401254 : 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 57796601 : gred_rfrac_simple(GEN n, GEN d)
260 : {
261 : GEN _1n, _1d, c, cn, cd, z;
262 57796601 : long dd = degpol(d);
263 :
264 57796601 : 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 57792533 : if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
272 57792533 : _1n = Rg_get_1(n);
273 57792533 : _1d = Rg_get_1(d);
274 57792533 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
275 57792533 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
276 57792526 : cd = content(d);
277 59662476 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
278 57792526 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
279 57792526 : if (!gequal1(cd)) {
280 6586874 : d = RgX_Rg_div(d,cd);
281 6586874 : if (!gequal1(cn))
282 : {
283 1313705 : 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 1313635 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
289 1313635 : c = gdiv(cn,cd);
290 : }
291 : }
292 : else
293 5273169 : c = ginv(cd);
294 : } else {
295 51205652 : 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 47864257 : GEN y = cgetg(3,t_RFRAC);
306 47864257 : gel(y,1) = gcopy(n);
307 47864257 : gel(y,2) = RgX_copy(d); return y;
308 : }
309 : }
310 :
311 9927233 : if (typ(c) == t_POL)
312 : {
313 912609 : z = c;
314 953335 : do { z = content(z); } while (typ(z) == t_POL);
315 912609 : cd = denom_i(z);
316 912609 : cn = gmul(c, cd);
317 : }
318 : else
319 : {
320 9014624 : cn = numer_i(c);
321 9014624 : cd = denom_i(c);
322 : }
323 9927233 : z = cgetg(3,t_RFRAC);
324 9927233 : gel(z,1) = gmul(n, cn);
325 9927233 : 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 9927233 : if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
328 9927233 : 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 44825298 : gred_rfrac2(GEN n, GEN d)
354 : {
355 : GEN y, z, _1n, _1d;
356 : long v, vd, vn;
357 :
358 44825298 : n = simplify_shallow(n);
359 44825298 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
360 37967947 : d = simplify_shallow(d);
361 37967947 : if (typ(d) != t_POL) return gdiv(n,d);
362 36791051 : vd = varn(d);
363 36791051 : 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 16134848 : vn = varn(n);
370 16134848 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
371 15992897 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
372 15833124 : _1n = Rg_get_1(n);
373 15833124 : _1d = Rg_get_1(d);
374 15833124 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
375 15833124 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
376 :
377 : /* now n and d are t_POLs in the same variable */
378 15833124 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
379 15833124 : 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 3035189 : if (!isinexact(n) && !isinexact(d))
387 : {
388 3034965 : y = RgX_divrem(n, d, &z);
389 3034965 : 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 142933695 : Qdivii(GEN x, GEN y)
399 : {
400 142933695 : pari_sp av = avma;
401 : GEN r, q;
402 :
403 142933695 : if (lgefint(y) == 3)
404 : {
405 124638018 : q = Qdiviu(x, y[2]);
406 124636118 : if (signe(y) > 0) return q;
407 11026168 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
408 11026426 : return q;
409 : }
410 18295677 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
411 18296872 : if (equali1(x))
412 : {
413 5186686 : if (!signe(y)) pari_err_INV("gdiv",y);
414 5186581 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
415 : }
416 13110202 : q = dvmdii(x,y,&r);
417 13110100 : if (r == gen_0) return q; /* gen_0 intended */
418 8660866 : r = gcdii(y, r);
419 8660812 : if (lgefint(r) == 3)
420 : {
421 7758313 : ulong t = r[2];
422 7758313 : set_avma(av);
423 7758290 : if (t == 1) q = mkfraccopy(x,y);
424 : else
425 : {
426 2799012 : q = cgetg(3,t_FRAC);
427 2799181 : gel(q,1) = diviuexact(x,t);
428 2799140 : gel(q,2) = diviuexact(y,t);
429 : }
430 : }
431 : else
432 : { /* rare: r and q left on stack for efficiency */
433 902499 : q = cgetg(3,t_FRAC);
434 902508 : gel(q,1) = diviiexact(x,r);
435 902518 : gel(q,2) = diviiexact(y,r);
436 : }
437 8660817 : normalize_frac(q); return q;
438 : }
439 :
440 : /* x t_INT, return x/y in reduced form */
441 : GEN
442 144127557 : Qdiviu(GEN x, ulong y)
443 : {
444 144127557 : pari_sp av = avma;
445 : ulong r, t;
446 : GEN q;
447 :
448 144127557 : if (y == 1) return icopy(x);
449 122219568 : if (!y) pari_err_INV("gdiv",gen_0);
450 122219568 : if (equali1(x)) retmkfrac(gen_1, utoipos(y));
451 115007962 : q = absdiviu_rem(x,y,&r);
452 115003608 : if (!r)
453 : {
454 70089469 : if (signe(x) < 0) togglesign(q);
455 70089385 : return q;
456 : }
457 44914139 : t = ugcd(y, r); set_avma(av);
458 44920316 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
459 17969697 : retmkfrac(diviuexact(x,t), utoipos(y / t));
460 : }
461 :
462 : /* x t_INT, return x/y in reduced form */
463 : GEN
464 1591257 : Qdivis(GEN x, long y)
465 : {
466 1591257 : pari_sp av = avma;
467 : ulong r, t;
468 : long s;
469 : GEN q;
470 :
471 1591257 : if (y > 0) return Qdiviu(x, y);
472 180383 : if (!y) pari_err_INV("gdiv",gen_0);
473 180383 : s = signe(x);
474 180383 : if (!s) return gen_0;
475 134064 : if (y < 0) { y = -y; s = -s; }
476 134064 : if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
477 133784 : if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
478 130928 : q = absdiviu_rem(x,y,&r);
479 130928 : if (!r)
480 : {
481 55860 : if (s < 0) togglesign(q);
482 55860 : return q;
483 : }
484 75068 : t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
485 75068 : if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
486 75068 : gel(q,1) = x; setsigne(x, s);
487 75068 : gel(q,2) = utoipos(y); return q;
488 : }
489 :
490 : /*******************************************************************/
491 : /* */
492 : /* CONJUGATION */
493 : /* */
494 : /*******************************************************************/
495 : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
496 : static GEN
497 18403 : quad_polmod_conj(GEN x, GEN y)
498 : {
499 : GEN z, u, v, a, b;
500 18403 : if (typ(x) != t_POL) return gcopy(x);
501 18403 : if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
502 18403 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
503 18403 : b = gel(y,3); v = gel(x,2);
504 18403 : z = cgetg(4, t_POL); z[1] = x[1];
505 18403 : gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
506 18403 : gel(z,3) = gneg(u); return z;
507 : }
508 : static GEN
509 18403 : quad_polmod_norm(GEN x, GEN y)
510 : {
511 : GEN z, u, v, a, b, c;
512 18403 : if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
513 18403 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
514 18403 : b = gel(y,3); v = gel(x,2);
515 18403 : c = gel(y,2);
516 18403 : z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
517 18403 : if (!gequal1(a)) z = gdiv(z, a);
518 18403 : return gadd(z, gsqr(v));
519 : }
520 :
521 : GEN
522 31554284 : conj_i(GEN x)
523 : {
524 31554284 : switch(typ(x))
525 : {
526 7461680 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
527 7461680 : return x;
528 :
529 23928484 : case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
530 :
531 945 : case t_QUAD:
532 : {
533 945 : GEN y = cgetg(4,t_QUAD);
534 945 : gel(y,1) = gel(x,1);
535 945 : gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
536 945 : : gadd(gel(x,2), gel(x,3));
537 945 : gel(y,3) = gneg(gel(x,3)); return y;
538 : }
539 350 : case t_POL: pari_APPLY_pol_normalized(conj_i(gel(x,i)));
540 31647 : case t_SER: pari_APPLY_ser_normalized(conj_i(gel(x,i)));
541 :
542 153440 : case t_RFRAC:
543 : case t_VEC:
544 : case t_COL:
545 553567 : case t_MAT: pari_APPLY_same(conj_i(gel(x,i)));
546 :
547 0 : case t_POLMOD:
548 : {
549 0 : GEN X = gel(x,1);
550 0 : long d = degpol(X);
551 0 : if (d < 2) return x;
552 0 : if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
553 : }
554 : }
555 0 : pari_err_TYPE("gconj",x);
556 : return NULL; /* LCOV_EXCL_LINE */
557 : }
558 : GEN
559 390135 : gconj(GEN x)
560 390135 : { pari_sp av = avma; return gc_GEN(av, conj_i(x)); }
561 :
562 : GEN
563 84 : conjvec(GEN x,long prec)
564 : {
565 : long lx, s, i;
566 : GEN z;
567 :
568 84 : switch(typ(x))
569 : {
570 0 : case t_INT: case t_INTMOD: case t_FRAC:
571 0 : return mkcolcopy(x);
572 :
573 0 : case t_COMPLEX: case t_QUAD:
574 0 : z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
575 :
576 28 : case t_FFELT:
577 28 : return FF_conjvec(x);
578 :
579 0 : case t_VEC: case t_COL:
580 0 : lx = lg(x); z = cgetg(lx,t_MAT);
581 0 : if (lx == 1) return z;
582 0 : gel(z,1) = conjvec(gel(x,1),prec);
583 0 : s = lgcols(z);
584 0 : for (i=2; i<lx; i++)
585 : {
586 0 : gel(z,i) = conjvec(gel(x,i),prec);
587 0 : if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
588 : }
589 0 : break;
590 :
591 56 : case t_POLMOD: {
592 56 : GEN T = gel(x,1), r;
593 : pari_sp av;
594 :
595 56 : lx = lg(T);
596 56 : if (lx <= 3) return cgetg(1,t_COL);
597 56 : x = gel(x,2);
598 238 : for (i=2; i<lx; i++)
599 : {
600 189 : GEN c = gel(T,i);
601 189 : switch(typ(c)) {
602 7 : case t_INTMOD: {
603 7 : GEN p = gel(c,1);
604 : pari_sp av;
605 7 : if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
606 7 : av = avma;
607 7 : T = RgX_to_FpX(T,p);
608 7 : x = RgX_to_FpX(x, p);
609 7 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
610 7 : z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
611 7 : return gc_upto(av, z);
612 : }
613 182 : case t_INT:
614 182 : case t_FRAC: break;
615 0 : default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
616 : }
617 : }
618 49 : if (typ(x) != t_POL)
619 : {
620 0 : if (!is_rational_t(typ(x)))
621 0 : pari_err_TYPE("conjvec [not a rational t_POL]",x);
622 0 : retconst_col(lx-3, gcopy(x));
623 : }
624 49 : RgX_check_QX(x,"conjvec");
625 49 : av = avma;
626 49 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
627 49 : r = cleanroots(T,prec);
628 49 : z = cgetg(lx-2,t_COL);
629 182 : for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
630 49 : return gc_upto(av, z);
631 : }
632 :
633 0 : default:
634 0 : pari_err_TYPE("conjvec",x);
635 : return NULL; /* LCOV_EXCL_LINE */
636 : }
637 0 : return z;
638 : }
639 :
640 : /********************************************************************/
641 : /** **/
642 : /** ADDITION **/
643 : /** **/
644 : /********************************************************************/
645 : static GEN
646 24592997 : mkpadic_mod(GEN u, GEN p, GEN pd, long e, long d)
647 24592997 : { retmkpadic(modii(u, pd), icopy(p), icopy(pd), e, d); }
648 :
649 : /* x, y compatible PADIC, op = add or sub */
650 : static GEN
651 19972825 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
652 : {
653 19972825 : pari_sp av = avma;
654 : long d, e, r, rx, ry;
655 19972825 : GEN u, z, mod, pdx, pdy, ux, uy, p = padic_p(x);
656 : int swap;
657 :
658 19972825 : e = valp(x);
659 19972825 : r = valp(y); d = r-e;
660 19972825 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
661 19972825 : pdx = padic_pd(x); ux = padic_u(x);
662 19972825 : pdy = padic_pd(y); uy = padic_u(y);
663 19972825 : rx = precp(x);
664 19972825 : ry = precp(y);
665 19972825 : if (d) /* v(x) < v(y) */
666 : {
667 10703707 : r = d+ry; z = powiu(p,d);
668 10703830 : if (r < rx) mod = mulii(z,pdy); else { r = rx; mod = pdx; }
669 10703834 : z = mulii(z, uy);
670 10703793 : u = swap? op(z, ux): op(ux, z);
671 : }
672 : else
673 : {
674 : long c;
675 9269118 : if (ry < rx) { r=ry; mod = pdy; } else { r=rx; mod = pdx; }
676 9269118 : u = swap? op(uy, ux): op(ux, uy);
677 9269615 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
678 : {
679 138287 : set_avma(av); return zeropadic(p, e+r);
680 : }
681 9131000 : if (c)
682 : {
683 3351555 : mod = diviiexact(mod, powiu(p,c));
684 3351553 : r -= c;
685 3351553 : e += c;
686 : }
687 : }
688 19834772 : return gc_upto(av, mkpadic_mod(u, p, mod, e, r));
689 : }
690 : /* Rg_to_Fp(t_FRAC) without GC */
691 : static GEN
692 226838 : Q_to_Fp(GEN x, GEN p)
693 226838 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
694 : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
695 : static GEN
696 4762997 : addQp(GEN x, GEN y)
697 : {
698 4762997 : pari_sp av = avma;
699 4762997 : long d, r, e, vy = valp(y), py = precp(y);
700 4762997 : GEN q, mod, u, p = padic_p(y);
701 :
702 4762997 : e = Q_pvalrem(x, p, &x);
703 4762979 : d = vy - e; r = d + py;
704 4762979 : if (r <= 0) { set_avma(av); return gcopy(y); }
705 4761070 : mod = padic_pd(y);
706 4761070 : u = padic_u(y);
707 4761070 : if (d > 0)
708 : {
709 1376553 : q = powiu(p,d);
710 1376560 : mod = mulii(mod, q);
711 1376561 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
712 1376561 : u = addii(x, mulii(u, q));
713 : }
714 3384517 : else if (d < 0)
715 : {
716 405332 : q = powiu(p,-d);
717 405332 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
718 405332 : u = addii(u, mulii(x, q));
719 405332 : r = py; e = vy;
720 : }
721 : else
722 : {
723 : long c;
724 2979185 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
725 2979185 : u = addii(u, x);
726 2979181 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
727 : {
728 1204 : set_avma(av); return zeropadic(p,e+r);
729 : }
730 2977965 : if (c)
731 : {
732 970644 : mod = diviiexact(mod, powiu(p,c));
733 970644 : r -= c;
734 970644 : e += c;
735 : }
736 : }
737 4759860 : return gc_upto(av, mkpadic_mod(u, p, mod, e, r));
738 : }
739 :
740 : /* Mod(x,X) + Mod(y,X) */
741 : #define addsub_polmod_same addsub_polmod_scal
742 : /* Mod(x,X) +/- Mod(y,Y) */
743 : static GEN
744 7203 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
745 : {
746 7203 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
747 7203 : GEN z = cgetg(3,t_POLMOD);
748 7203 : long vx = varn(X), vy = varn(Y);
749 7203 : if (vx==vy) {
750 : pari_sp av;
751 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
752 14 : warn_coercion(X,Y,gel(z,1));
753 14 : gel(z,2) = gc_upto(av, gmod(op(x, y), gel(z,1))); return z;
754 : }
755 7189 : if (varncmp(vx, vy) < 0)
756 7189 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
757 : else
758 0 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
759 7189 : gel(z,2) = op(x, y); return z;
760 : }
761 : /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
762 : static GEN
763 13455136 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
764 13455136 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
765 :
766 : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
767 : static GEN
768 408206 : add_ser_scal(GEN y, GEN x)
769 : {
770 : long i, v, ly, vy;
771 : GEN z;
772 :
773 408206 : if (isrationalzero(x)) return gcopy(y);
774 381921 : ly = lg(y);
775 381921 : v = valser(y);
776 381921 : if (v < 3-ly) return gcopy(y);
777 : /* v + ly >= 3 */
778 381669 : if (v < 0)
779 : {
780 1162 : z = cgetg(ly,t_SER); z[1] = y[1];
781 3276 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
782 1162 : gel(z,i) = gadd(x,gel(y,i)); i++;
783 3157 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
784 1162 : return normalizeser(z);
785 : }
786 380507 : vy = varn(y);
787 380507 : if (v > 0)
788 : {
789 19355 : if (ser_isexactzero(y))
790 9394 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
791 9961 : y -= v; ly += v;
792 9961 : z = cgetg(ly,t_SER);
793 9961 : x = gcopy(x);
794 20559 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
795 : }
796 361152 : else if (ser_isexactzero(y)) return gcopy(y);
797 : else
798 : { /* v = 0, ly >= 3 */
799 361145 : z = cgetg(ly,t_SER);
800 361145 : x = gadd(x, gel(y,2));
801 361145 : i = 3;
802 : }
803 1599007 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
804 371106 : gel(z,2) = x;
805 371106 : z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
806 371106 : return gequal0(x)? normalizeser(z): z;
807 : }
808 : static long
809 7162062 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
810 : /* x,y t_SER in the same variable: x+y */
811 : static GEN
812 3581430 : ser_add(GEN x, GEN y)
813 : {
814 3581430 : long i, lx,ly, n = valser(y) - valser(x);
815 : GEN z;
816 3581430 : if (n < 0) { n = -n; swap(x,y); }
817 : /* valser(x) <= valser(y) */
818 3581430 : lx = _serprec(x);
819 3581430 : if (lx == 2) /* don't lose type information */
820 : {
821 798 : z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
822 798 : setvalser(z, valser(x)); return z;
823 : }
824 3580632 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
825 3580632 : if (n)
826 : {
827 107803 : if (n+2 > lx) return gcopy(x);
828 107103 : z = cgetg(ly,t_SER);
829 819228 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
830 506771 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
831 : } else {
832 3472829 : z = cgetg(ly,t_SER);
833 17589564 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
834 : }
835 3579932 : z[1] = x[1]; return normalizeser(z);
836 : }
837 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
838 : static GEN
839 8811954 : add_rfrac_scal(GEN y, GEN x)
840 : {
841 : pari_sp av;
842 : GEN n;
843 :
844 8811954 : if (isintzero(x)) return gcopy(y); /* frequent special case */
845 5138189 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
846 5138189 : return gc_upto(av, gred_rfrac_simple(n, gel(y,2)));
847 : }
848 :
849 : /* x "scalar", ty != t_MAT and nonscalar */
850 : static GEN
851 45452063 : add_scal(GEN y, GEN x, long ty)
852 : {
853 45452063 : switch(ty)
854 : {
855 40529624 : case t_POL: return RgX_Rg_add(y, x);
856 408178 : case t_SER: return add_ser_scal(y, x);
857 4468855 : case t_RFRAC: return add_rfrac_scal(y, x);
858 0 : case t_COL: return RgC_Rg_add(y, x);
859 45399 : case t_VEC:
860 45399 : if (isintzero(x)) return gcopy(y);
861 175 : break;
862 : }
863 182 : pari_err_TYPE2("+",x,y);
864 : return NULL; /* LCOV_EXCL_LINE */
865 : }
866 :
867 : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
868 : static GEN
869 15041289 : setfrac(GEN z, GEN a, GEN b)
870 : {
871 15041289 : gel(z,1) = icopy_avma(a, (pari_sp)z);
872 15041283 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
873 15041347 : set_avma((pari_sp)gel(z,2)); return z;
874 : }
875 : /* z <- a / (b*Q), (Q,a) = 1 */
876 : static GEN
877 14226515 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
878 : {
879 14226515 : GEN q = Qdivii(a, b); /* != 0 */
880 14226584 : if (typ(q) == t_INT)
881 : {
882 1960475 : gel(z,1) = gc_INT((pari_sp)Q, q);
883 1960475 : gel(z,2) = Q; return z;
884 : }
885 12266109 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
886 : }
887 : static GEN
888 42657411 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
889 : {
890 42657411 : GEN x1 = gel(x,1), x2 = gel(x,2);
891 42657411 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
892 42657411 : int s = cmpii(x2, y2);
893 :
894 : /* common denominator: (x1 op y1) / x2 */
895 42657404 : if (!s)
896 : {
897 17110114 : pari_sp av = avma;
898 17110114 : return gc_upto(av, Qdivii(op(x1, y1), x2));
899 : }
900 25547290 : z = cgetg(3, t_FRAC);
901 25547307 : if (s < 0)
902 : {
903 17706225 : Q = dvmdii(y2, x2, &r);
904 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
905 17706195 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
906 9752940 : delta = gcdii(x2,r);
907 : }
908 : else
909 : {
910 7841082 : Q = dvmdii(x2, y2, &r);
911 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
912 7841096 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
913 1567769 : delta = gcdii(y2,r);
914 : }
915 : /* delta = gcd(x2,y2) */
916 11320733 : if (equali1(delta))
917 : { /* numerator is nonzero */
918 8545477 : gel(z,1) = gc_INT((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
919 8545477 : gel(z,2) = mulii(x2,y2); return z;
920 : }
921 2775254 : x2 = diviiexact(x2,delta);
922 2775258 : y2 = diviiexact(y2,delta);
923 2775256 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
924 2775257 : q = dvmdii(n, delta, &r);
925 2775253 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
926 2521274 : r = gcdii(delta, r);
927 2521281 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
928 2521281 : return setfrac(z, n, mulii(mulii(x2, y2), delta));
929 : }
930 :
931 : /* assume x2, y2 are t_POLs in the same variable */
932 : static GEN
933 3018271 : add_rfrac(GEN x, GEN y)
934 : {
935 3018271 : pari_sp av = avma;
936 3018271 : GEN x1 = gel(x,1), x2 = gel(x,2);
937 3018271 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
938 :
939 3018271 : delta = RgX_gcd(x2,y2);
940 3018271 : if (!degpol(delta))
941 : {
942 644 : n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
943 644 : d = RgX_mul(x2, y2);
944 644 : return gc_upto(av, gred_rfrac_simple(n, d));
945 : }
946 3017627 : x2 = RgX_div(x2,delta);
947 3017627 : y2 = RgX_div(y2,delta);
948 3017627 : n = gadd(gmul(x1,y2), gmul(y1,x2));
949 3017627 : if (!signe(n))
950 : {
951 721418 : n = simplify_shallow(n);
952 721418 : if (isexactzero(n))
953 : {
954 721411 : if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
955 7 : return gc_upto(av, scalarpol(n, varn(x2)));
956 : }
957 7 : return gc_GEN(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
958 : }
959 2296209 : if (degpol(n) == 0)
960 1150575 : return gc_upto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
961 1145634 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
962 1145634 : if (isexactzero(r))
963 : {
964 : GEN z;
965 227982 : d = RgX_mul(x2, y2);
966 : /* "constant" denominator ? */
967 227982 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
968 227982 : return gc_upto(av, z);
969 : }
970 917652 : r = RgX_gcd(delta, r);
971 917652 : if (degpol(r))
972 : {
973 160642 : n = RgX_div(n, r);
974 160642 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
975 : }
976 : else
977 757010 : d = RgX_mul(gel(x,2), y2);
978 917652 : return gc_upto(av, gred_rfrac_simple(n, d));
979 : }
980 :
981 : GEN
982 6111045660 : gadd(GEN x, GEN y)
983 : {
984 6111045660 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
985 : pari_sp av;
986 : GEN z, p1;
987 :
988 6111045660 : if (tx == ty) switch(tx) /* shortcut to generic case */
989 : {
990 3004075355 : case t_INT: return addii(x,y);
991 1654047541 : case t_REAL: return addrr(x,y);
992 1176509 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
993 1176509 : z = cgetg(3,t_INTMOD);
994 1176509 : if (X==Y || equalii(X,Y))
995 1176495 : return add_intmod_same(z, X, gel(x,2), gel(y,2));
996 14 : gel(z,1) = gcdii(X,Y);
997 14 : warn_coercion(X,Y,gel(z,1));
998 14 : av = avma; p1 = addii(gel(x,2),gel(y,2));
999 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
1000 : }
1001 36584486 : case t_FRAC: return addsub_frac(x,y,addii);
1002 358650400 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1003 358156940 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1004 358836785 : if (isintzero(gel(z,2)))
1005 : {
1006 521178 : set_avma((pari_sp)(z+3));
1007 521178 : return gadd(gel(x,1),gel(y,1));
1008 : }
1009 358213356 : gel(z,1) = gadd(gel(x,1),gel(y,1));
1010 358223177 : return z;
1011 14530642 : case t_PADIC:
1012 14530642 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1013 14530248 : return addsub_pp(x,y, addii);
1014 672 : case t_QUAD: z = cgetg(4,t_QUAD);
1015 672 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1016 672 : gel(z,1) = ZX_copy(gel(x,1));
1017 672 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1018 672 : gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
1019 10151518 : case t_POLMOD:
1020 10151518 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1021 10144350 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
1022 7168 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
1023 13133909 : case t_FFELT: return FF_add(x,y);
1024 140366255 : case t_POL:
1025 140366255 : vx = varn(x);
1026 140366255 : vy = varn(y);
1027 140366255 : if (vx != vy) {
1028 4863731 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1029 1801515 : else return RgX_Rg_add(y, x);
1030 : }
1031 135502524 : return RgX_add(x, y);
1032 3576068 : case t_SER:
1033 3576068 : vx = varn(x);
1034 3576068 : vy = varn(y);
1035 3576068 : if (vx != vy) {
1036 28 : if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1037 21 : else return add_ser_scal(y, x);
1038 : }
1039 3576040 : return ser_add(x, y);
1040 4310172 : case t_RFRAC:
1041 4310172 : vx = varn(gel(x,2));
1042 4310172 : vy = varn(gel(y,2));
1043 4310172 : if (vx != vy) {
1044 1291901 : if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1045 538397 : else return add_rfrac_scal(y, x);
1046 : }
1047 3018271 : return add_rfrac(x,y);
1048 2115988 : case t_VEC:
1049 2115988 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1050 2115988 : return RgV_add(x,y);
1051 1113399 : case t_COL:
1052 1113399 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1053 1113399 : return RgC_add(x,y);
1054 671412 : case t_MAT:
1055 671412 : lx = lg(x);
1056 671412 : if (lg(y) != lx) pari_err_OP("+",x,y);
1057 671412 : if (lx == 1) return cgetg(1, t_MAT);
1058 671412 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1059 671405 : return RgM_add(x,y);
1060 :
1061 0 : default: pari_err_TYPE2("+",x,y);
1062 : }
1063 : /* tx != ty */
1064 870813249 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1065 :
1066 870813249 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1067 : {
1068 753582257 : case t_INT:
1069 : switch(ty)
1070 : {
1071 440282437 : case t_REAL: return addir(x,y);
1072 2138731 : case t_INTMOD:
1073 2138731 : z = cgetg(3, t_INTMOD);
1074 2138731 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1075 51809891 : case t_FRAC: z = cgetg(3,t_FRAC);
1076 51809887 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1077 51809825 : gel(z,2) = icopy(gel(y,2)); return z;
1078 252433672 : case t_COMPLEX: return addRc(x, y);
1079 5597511 : case t_PADIC:
1080 5597511 : if (!signe(x)) return gcopy(y);
1081 4533227 : return addQp(x,y);
1082 1113 : case t_QUAD: return addRq(x, y);
1083 1365044 : case t_FFELT: return FF_Z_add(y,x);
1084 : }
1085 :
1086 : case t_REAL:
1087 58321966 : switch(ty)
1088 : {
1089 14388666 : case t_FRAC:
1090 14388666 : if (!signe(gel(y,1))) return rcopy(x);
1091 14388666 : if (!signe(x))
1092 : {
1093 12080 : lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1094 12080 : return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1095 : }
1096 14376586 : av = avma; z = addir(gel(y,1), mulir(gel(y,2),x));
1097 14374983 : return gc_leaf(av, divri(z,gel(y,2)));
1098 43933230 : case t_COMPLEX: return addRc(x, y);
1099 70 : case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, realprec(x));
1100 :
1101 0 : default: pari_err_TYPE2("+",x,y);
1102 : }
1103 :
1104 17647 : case t_INTMOD:
1105 : switch(ty)
1106 : {
1107 17507 : case t_FRAC: { GEN X = gel(x,1);
1108 17507 : z = cgetg(3, t_INTMOD);
1109 17507 : p1 = Fp_div(gel(y,1), gel(y,2), X);
1110 17507 : return add_intmod_same(z, X, p1, gel(x,2));
1111 : }
1112 14 : case t_FFELT:
1113 14 : if (!equalii(gel(x,1),FF_p_i(y)))
1114 0 : pari_err_OP("+",x,y);
1115 14 : return FF_Z_add(y,gel(x,2));
1116 91 : case t_COMPLEX: return addRc(x, y);
1117 0 : case t_PADIC: { GEN X = gel(x,1);
1118 0 : z = cgetg(3, t_INTMOD);
1119 0 : return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1120 : }
1121 35 : case t_QUAD: return addRq(x, y);
1122 : }
1123 :
1124 : case t_FRAC:
1125 : switch (ty)
1126 : {
1127 10317828 : case t_COMPLEX: return addRc(x, y);
1128 229771 : case t_PADIC:
1129 229771 : if (!signe(gel(x,1))) return gcopy(y);
1130 229771 : return addQp(x,y);
1131 133 : case t_QUAD: return addRq(x, y);
1132 1337 : case t_FFELT: return FF_Q_add(y, x);
1133 : }
1134 :
1135 : case t_FFELT:
1136 0 : pari_err_TYPE2("+",x,y);
1137 :
1138 35 : case t_COMPLEX:
1139 : switch(ty)
1140 : {
1141 28 : case t_PADIC:
1142 28 : return Zp_nosquare_m1(padic_p(y))? addRc(y, x): addTp(x, y);
1143 7 : case t_QUAD:
1144 7 : lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1145 7 : return gequal0(y)? gcopy(x): addqf(y, x, lx);
1146 : }
1147 :
1148 : case t_PADIC: /* ty == t_QUAD */
1149 6 : return (kro_quad(y, padic_p(x)) == -1)? addRq(x, y): addTp(y, x);
1150 : }
1151 : /* tx < ty, !is_const_t(y) */
1152 51143594 : switch(ty)
1153 : {
1154 7693 : case t_MAT:
1155 7693 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1156 7693 : if (isrationalzero(x)) return gcopy(y);
1157 7616 : return RgM_Rg_add(y, x);
1158 186395 : case t_COL:
1159 186395 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1160 186395 : return RgC_Rg_add(y, x);
1161 2409753 : case t_POLMOD: /* is_const_t(tx) in this case */
1162 2409753 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1163 : }
1164 48539753 : if (is_scalar_t(tx)) {
1165 45459893 : if (tx == t_POLMOD)
1166 : {
1167 117810 : vx = varn(gel(x,1));
1168 117810 : vy = gvar(y);
1169 117810 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1170 : else
1171 87759 : if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1172 30590 : return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1173 : }
1174 45342083 : return add_scal(y, x, ty);
1175 : }
1176 : /* x and y are not scalars, ty != t_MAT */
1177 3079870 : vx = gvar(x);
1178 3079872 : vy = gvar(y);
1179 3079872 : if (vx != vy) { /* x or y is treated as a scalar */
1180 22759 : if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1181 32417 : return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1182 32417 : : add_scal(y, x, ty);
1183 : }
1184 : /* vx = vy */
1185 3057113 : switch(tx)
1186 : {
1187 3056616 : case t_POL:
1188 : switch (ty)
1189 : {
1190 5418 : case t_SER:
1191 5418 : if (lg(x) == 2) return gcopy(y);
1192 5397 : i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1193 5397 : i = lg(y) + valser(y) - i;
1194 5397 : if (i < 3) return gcopy(y);
1195 5390 : p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1196 5390 : settyp(p1, t_VECSMALL); /* p1 left on stack */
1197 5390 : return y;
1198 :
1199 3051198 : case t_RFRAC: return add_rfrac_scal(y, x);
1200 : }
1201 0 : break;
1202 :
1203 497 : case t_SER:
1204 497 : if (ty == t_RFRAC)
1205 : {
1206 : long vn, vd;
1207 497 : av = avma;
1208 497 : vn = gval(gel(y,1), vy);
1209 497 : vd = RgX_valrem_inexact(gel(y,2), NULL);
1210 497 : if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
1211 :
1212 497 : l = lg(x) + valser(x) - (vn - vd);
1213 497 : if (l < 3) { set_avma(av); return gcopy(x); }
1214 497 : return gc_upto(av, gadd(x, rfrac_to_ser_i(y, l)));
1215 : }
1216 0 : break;
1217 : }
1218 0 : pari_err_TYPE2("+",x,y);
1219 : return NULL; /* LCOV_EXCL_LINE */
1220 : }
1221 :
1222 : GEN
1223 270736815 : gaddsg(long x, GEN y)
1224 : {
1225 270736815 : long ty = typ(y);
1226 : GEN z;
1227 :
1228 270736815 : switch(ty)
1229 : {
1230 124844288 : case t_INT: return addsi(x,y);
1231 119268571 : case t_REAL: return addsr(x,y);
1232 12061 : case t_INTMOD:
1233 12061 : z = cgetg(3, t_INTMOD);
1234 12061 : return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1235 15248671 : case t_FRAC: z = cgetg(3,t_FRAC);
1236 15248671 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1237 15248671 : gel(z,2) = icopy(gel(y,2)); return z;
1238 8225328 : case t_COMPLEX:
1239 8225328 : z = cgetg(3, t_COMPLEX);
1240 8225328 : gel(z,1) = gaddsg(x, gel(y,1));
1241 8225328 : gel(z,2) = gcopy(gel(y,2)); return z;
1242 :
1243 3137896 : default: return gadd(stoi(x), y);
1244 : }
1245 : }
1246 :
1247 : GEN
1248 3238194 : gsubsg(long x, GEN y)
1249 : {
1250 : GEN z, a, b;
1251 : pari_sp av;
1252 :
1253 3238194 : switch(typ(y))
1254 : {
1255 276787 : case t_INT: return subsi(x,y);
1256 1234941 : case t_REAL: return subsr(x,y);
1257 56 : case t_INTMOD:
1258 56 : z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1259 56 : return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1260 732506 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1261 732506 : gel(z,1) = gc_INT((pari_sp)z, subii(mulis(b,x), a));
1262 732506 : gel(z,2) = icopy(gel(y,2)); return z;
1263 956790 : case t_COMPLEX:
1264 956790 : z = cgetg(3, t_COMPLEX);
1265 956790 : gel(z,1) = gsubsg(x, gel(y,1));
1266 956790 : gel(z,2) = gneg(gel(y,2)); return z;
1267 : }
1268 37114 : av = avma;
1269 37114 : return gc_upto(av, gadd(stoi(x), gneg_i(y)));
1270 : }
1271 :
1272 : /********************************************************************/
1273 : /** **/
1274 : /** SUBTRACTION **/
1275 : /** **/
1276 : /********************************************************************/
1277 :
1278 : GEN
1279 2877140165 : gsub(GEN x, GEN y)
1280 : {
1281 2877140165 : long tx = typ(x), ty = typ(y);
1282 : pari_sp av;
1283 : GEN z;
1284 2877140165 : if (tx == ty) switch(tx) /* shortcut to generic case */
1285 : {
1286 2086529830 : case t_INT: return subii(x,y);
1287 585578215 : case t_REAL: return subrr(x,y);
1288 1027943 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1289 1027943 : z = cgetg(3,t_INTMOD);
1290 1027943 : if (X==Y || equalii(X,Y))
1291 1027929 : return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1292 14 : gel(z,1) = gcdii(X,Y);
1293 14 : warn_coercion(X,Y,gel(z,1));
1294 14 : av = avma; p1 = subii(gel(x,2),gel(y,2));
1295 14 : gel(z,2) = gc_INT(av, modii(p1, gel(z,1))); return z;
1296 : }
1297 6072939 : case t_FRAC: return addsub_frac(x,y, subii);
1298 102268327 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1299 102310396 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1300 102209375 : if (isintzero(gel(z,2)))
1301 : {
1302 21322 : set_avma((pari_sp)(z+3));
1303 21322 : return gsub(gel(x,1),gel(y,1));
1304 : }
1305 102193932 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1306 102174572 : return z;
1307 5442665 : case t_PADIC:
1308 5442665 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1309 5442665 : return addsub_pp(x,y, subii);
1310 168 : case t_QUAD: z = cgetg(4,t_QUAD);
1311 168 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1312 168 : gel(z,1) = ZX_copy(gel(x,1));
1313 168 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1314 168 : gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1315 870478 : case t_POLMOD:
1316 870478 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1317 870443 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1318 35 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1319 240348 : case t_FFELT: return FF_sub(x,y);
1320 20214549 : case t_POL: {
1321 20214549 : long vx = varn(x);
1322 20214549 : long vy = varn(y);
1323 20214549 : if (vx != vy) {
1324 1150579 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1325 85352 : else return Rg_RgX_sub(x, y);
1326 : }
1327 19063970 : return RgX_sub(x, y);
1328 : }
1329 299410 : case t_VEC:
1330 299410 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1331 299410 : return RgV_sub(x,y);
1332 3581714 : case t_COL:
1333 3581714 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1334 3581714 : return RgC_sub(x,y);
1335 254035 : case t_MAT: {
1336 254035 : long lx = lg(x);
1337 254035 : if (lg(y) != lx) pari_err_OP("+",x,y);
1338 254037 : if (lx == 1) return cgetg(1, t_MAT);
1339 173945 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1340 173943 : return RgM_sub(x,y);
1341 : }
1342 2448387 : case t_RFRAC: case t_SER: break;
1343 :
1344 0 : default: pari_err_TYPE2("+",x,y);
1345 : }
1346 66433058 : av = avma;
1347 66433058 : return gc_upto(av, gadd(x,gneg_i(y)));
1348 : }
1349 :
1350 : /********************************************************************/
1351 : /** **/
1352 : /** MULTIPLICATION **/
1353 : /** **/
1354 : /********************************************************************/
1355 : static GEN
1356 316386 : mul_ser_scal(GEN x, GEN t)
1357 : {
1358 316386 : if (isexactzero(t)) return gmul(Rg_get_0(x), t);
1359 313138 : if (isint1(t)) return gcopy(x);
1360 262514 : if (ser_isexactzero(x))
1361 : {
1362 378 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gmul(gel(x,2), t), varn(x), 1);
1363 378 : setvalser(z, valser(x)); return z;
1364 : }
1365 3442782 : pari_APPLY_ser(gmul(gel(x,i), t));
1366 : }
1367 : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1368 : * [n/d a valid RFRAC] */
1369 : static GEN
1370 10454828 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1371 : {
1372 10454828 : pari_sp av = avma;
1373 : GEN z;
1374 :
1375 10454828 : switch(typ(x))
1376 : {
1377 21 : case t_PADIC:
1378 21 : n = gmul(n, x);
1379 21 : d = gcvtop(d, padic_p(x), signe(padic_u(x))? precp(x): 1);
1380 21 : return gc_upto(av, gdiv(n,d));
1381 :
1382 630 : case t_INTMOD: case t_POLMOD:
1383 630 : n = gmul(n, x);
1384 630 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1385 630 : return gc_upto(av, gdiv(n,d));
1386 : }
1387 10454177 : z = gred_rfrac2(x, d);
1388 10454177 : n = simplify_shallow(n);
1389 10454177 : if (typ(z) == t_RFRAC)
1390 : {
1391 7929477 : n = gmul(gel(z,1), n);
1392 7929477 : d = gel(z,2);
1393 7929477 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1394 0 : z = RgX_Rg_div(n, d);
1395 : else
1396 7929477 : z = gred_rfrac_simple(n, d);
1397 : }
1398 : else
1399 2524700 : z = gmul(z, n);
1400 10454177 : return gc_upto(av, z);
1401 : }
1402 : static GEN
1403 157513740 : mul_scal(GEN y, GEN x, long ty)
1404 : {
1405 157513740 : switch(ty)
1406 : {
1407 148790778 : case t_POL:
1408 148790778 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1409 136725033 : return RgX_Rg_mul(y, x);
1410 308245 : case t_SER: return mul_ser_scal(y, x);
1411 8414709 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1412 14 : case t_QFB:
1413 14 : if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1414 : }
1415 12 : pari_err_TYPE2("*",x,y);
1416 : return NULL; /* LCOV_EXCL_LINE */
1417 : }
1418 : static GEN
1419 7545773 : mul_self_scal(GEN x, GEN y)
1420 634441579 : { pari_APPLY_same(gmul(y,gel(x,i))); }
1421 :
1422 : static GEN
1423 161432 : mul_gen_rfrac(GEN X, GEN Y)
1424 : {
1425 161432 : GEN y1 = gel(Y,1), y2 = gel(Y,2);
1426 161432 : long vx = gvar(X), vy = varn(y2);
1427 167193 : return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1428 5761 : gred_rfrac_simple(gmul(y1,X), y2);
1429 : }
1430 : /* (x1/x2) * (y1/y2) */
1431 : static GEN
1432 7910285 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1433 : {
1434 : GEN z, X, Y;
1435 7910285 : pari_sp av = avma;
1436 :
1437 7910285 : X = gred_rfrac2(x1, y2);
1438 7910285 : Y = gred_rfrac2(y1, x2);
1439 7910285 : if (typ(X) == t_RFRAC)
1440 : {
1441 6631189 : if (typ(Y) == t_RFRAC) {
1442 6565032 : x1 = gel(X,1);
1443 6565032 : x2 = gel(X,2);
1444 6565032 : y1 = gel(Y,1);
1445 6565032 : y2 = gel(Y,2);
1446 6565032 : z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1447 : } else
1448 66157 : z = mul_gen_rfrac(Y, X);
1449 : }
1450 1279096 : else if (typ(Y) == t_RFRAC)
1451 95275 : z = mul_gen_rfrac(X, Y);
1452 : else
1453 1183821 : z = gmul(X, Y);
1454 7910285 : return gc_upto(av, z);
1455 : }
1456 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1457 : static GEN
1458 262344 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1459 : {
1460 262344 : pari_sp av = avma;
1461 262344 : GEN X = gred_rfrac2(x1, y2);
1462 262344 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1463 : {
1464 255764 : x2 = RgX_mul(gel(X,2), x2);
1465 255764 : x1 = gel(X,1);
1466 : }
1467 : else
1468 6580 : x1 = X;
1469 262344 : return gc_upto(av, gred_rfrac_simple(x1, x2));
1470 : }
1471 :
1472 : /* Mod(y, Y) * x, assuming x scalar */
1473 : static GEN
1474 3238383 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1475 3238383 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1476 :
1477 : /* cf mulqq */
1478 : static GEN
1479 5860097 : quad_polmod_mul(GEN T, GEN x, GEN y)
1480 : {
1481 5860097 : GEN b = gel(T,3), c = gel(T,2);
1482 5860097 : GEN ux = gel(x,2), vx = gel(x,3), uy = gel(y,2), vy = gel(y,3);
1483 : GEN z, s, U, V, E, F;
1484 : pari_sp av, av2;
1485 5860097 : z = cgetg(4, t_POL); z[1] = x[1]; av = avma;
1486 5860097 : U = gmul(ux, uy);
1487 5860097 : V = gmul(vx, vy); s = gmul(c, V);
1488 5860097 : E = gmul(gadd(ux, vx), gadd(uy, vy));
1489 5860097 : if (typ(b) != t_INT) F = gadd(U, gmul(gaddgs(b, 1), V)); /* = U + (b+1) V */
1490 : else
1491 : { /* minor optimization */
1492 5860076 : if (!signe(b)) F = gadd(U, V); /* b = 0 */
1493 4284989 : else if (is_pm1(b)) /* b = 1 or -1 */
1494 4284331 : F = signe(b) < 0? U: gadd(U, gmul2n(V,1));
1495 : else /* |b| > 1 */
1496 658 : F = gadd(U, gmul(addis(b, 1), V));
1497 : }
1498 5860097 : av2 = avma;
1499 5860097 : gel(z,2) = gsub(U, s);
1500 5860097 : gel(z,3) = gsub(E, F);
1501 5860097 : gc_slice_unsafe(av, av2, z+2, 2);
1502 5860097 : return normalizepol_lg(z,4);
1503 : }
1504 : /* Mod(x,T) * Mod(y,T) */
1505 : static GEN
1506 8528815 : mul_polmod_same(GEN T, GEN x, GEN y)
1507 : {
1508 8528815 : GEN z = cgetg(3,t_POLMOD), a;
1509 8528815 : long v = varn(T), lx = lg(x), ly = lg(y);
1510 8528815 : gel(z,1) = RgX_copy(T);
1511 : /* x * y mod T optimised */
1512 8528815 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1513 7902082 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1514 1545858 : a = gmul(x, y);
1515 : else
1516 : {
1517 6982957 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1518 5854819 : a = quad_polmod_mul(T, x, y);
1519 : else
1520 1128138 : a = RgXQ_mul(x, y, gel(z,1));
1521 : }
1522 8528815 : gel(z,2) = a; return z;
1523 : }
1524 : static GEN
1525 419385 : sqr_polmod(GEN T, GEN x)
1526 : {
1527 419385 : GEN a, z = cgetg(3,t_POLMOD);
1528 419385 : gel(z,1) = RgX_copy(T);
1529 419385 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1530 28967 : a = gsqr(x);
1531 : else
1532 : {
1533 390418 : pari_sp av = avma;
1534 390418 : a = RgXQ_sqr(x, gel(z,1));
1535 390418 : a = gc_upto(av, a);
1536 : }
1537 419385 : gel(z,2) = a; return z;
1538 : }
1539 : /* Mod(x,X) * Mod(y,Y) */
1540 : static GEN
1541 2674945 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1542 : {
1543 2674945 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1544 2674945 : long vx = varn(X), vy = varn(Y);
1545 2674945 : GEN z = cgetg(3,t_POLMOD);
1546 :
1547 2674945 : if (vx==vy) {
1548 : pari_sp av;
1549 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
1550 14 : warn_coercion(X,Y,gel(z,1));
1551 14 : gel(z,2) = gc_upto(av, gmod(gmul(x, y), gel(z,1)));
1552 14 : return z;
1553 : }
1554 2674931 : if (varncmp(vx, vy) < 0)
1555 410662 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1556 : else
1557 2264269 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1558 2674931 : gel(z,2) = gmul(x, y); return z;
1559 : }
1560 :
1561 : #if 0 /* used by 3M only */
1562 : /* set z = x+y and return 1 if x,y have the same sign
1563 : * set z = x-y and return 0 otherwise */
1564 : static int
1565 : did_add(GEN x, GEN y, GEN *z)
1566 : {
1567 : long tx = typ(x), ty = typ(y);
1568 : if (tx == ty) switch(tx)
1569 : {
1570 : case t_INT: *z = addii(x,y); return 1;
1571 : case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1572 : case t_REAL:
1573 : if (signe(x) == -signe(y))
1574 : { *z = subrr(x,y); return 0; }
1575 : else
1576 : { *z = addrr(x,y); return 1; }
1577 : }
1578 : if (tx == t_REAL) switch(ty)
1579 : {
1580 : case t_INT:
1581 : if (signe(x) == -signe(y))
1582 : { *z = subri(x,y); return 0; }
1583 : else
1584 : { *z = addri(x,y); return 1; }
1585 : case t_FRAC:
1586 : if (signe(x) == -signe(gel(y,1)))
1587 : { *z = gsub(x,y); return 0; }
1588 : else
1589 : { *z = gadd(x,y); return 1; }
1590 : }
1591 : else if (ty == t_REAL) switch(tx)
1592 : {
1593 : case t_INT:
1594 : if (signe(x) == -signe(y))
1595 : { *z = subir(x,y); return 0; }
1596 : else
1597 : { *z = addir(x,y); return 1; }
1598 : case t_FRAC:
1599 : if (signe(gel(x,1)) == -signe(y))
1600 : { *z = gsub(x,y); return 0; }
1601 : else
1602 : { *z = gadd(x,y); return 1; }
1603 : }
1604 : *z = gadd(x,y); return 1;
1605 : }
1606 : #endif
1607 : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1608 : static GEN
1609 11668883 : mulcIR(GEN x, GEN y)
1610 : {
1611 11668883 : GEN z = cgetg(3,t_COMPLEX);
1612 11669064 : pari_sp av = avma;
1613 11669064 : gel(z,1) = gc_upto(av, gneg(gmul(y,gel(x,2))));
1614 11669254 : gel(z,2) = gmul(y, gel(x,1));
1615 11669114 : return z;
1616 :
1617 : }
1618 : /* x,y COMPLEX */
1619 : static GEN
1620 278790444 : mulcc(GEN x, GEN y)
1621 : {
1622 278790444 : GEN xr = gel(x,1), xi = gel(x,2);
1623 278790444 : GEN yr = gel(y,1), yi = gel(y,2);
1624 : GEN p1, p2, p3, p4, z;
1625 : pari_sp tetpil, av;
1626 :
1627 278790444 : if (isintzero(xr))
1628 : {
1629 15650303 : if (isintzero(yr)) {
1630 7129317 : av = avma;
1631 7129317 : return gc_upto(av, gneg(gmul(xi,yi)));
1632 : }
1633 8520819 : return mulcIR(y, xi);
1634 : }
1635 263130231 : if (isintzero(yr)) return mulcIR(x, yi);
1636 :
1637 259973960 : z = cgetg(3,t_COMPLEX); av = avma;
1638 : #if 0
1639 : /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1640 : * e.g. xr + xi if exponents differ */
1641 : if (did_add(xr, xi, &p3))
1642 : {
1643 : if (did_add(yr, yi, &p4)) {
1644 : /* R = xr*yr - xi*yi
1645 : * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1646 : p1 = gmul(xr,yr);
1647 : p2 = gmul(xi,yi); p2 = gneg(p2);
1648 : p3 = gmul(p3, p4);
1649 : p4 = gsub(p2, p1);
1650 : } else {
1651 : /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1652 : * I = xr*yi + xi*yr */
1653 : p1 = gmul(p3,p4);
1654 : p3 = gmul(xr,yi);
1655 : p4 = gmul(xi,yr);
1656 : p2 = gsub(p3, p4);
1657 : }
1658 : } else {
1659 : if (did_add(yr, yi, &p4)) {
1660 : /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1661 : * I = xr*yi +xi*yr */
1662 : p1 = gmul(p3,p4);
1663 : p3 = gmul(xr,yi);
1664 : p4 = gmul(xi,yr);
1665 : p2 = gsub(p4, p3);
1666 : } else {
1667 : /* R = xr*yr - xi*yi
1668 : * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1669 : p3 = gneg( gmul(p3, p4) );
1670 : p1 = gmul(xr,yr);
1671 : p2 = gmul(xi,yi);
1672 : p4 = gadd(p1, p2);
1673 :
1674 : p2 = gneg(p2);
1675 : }
1676 : }
1677 : tetpil = avma;
1678 : gel(z,1) = gadd(p1,p2);
1679 : gel(z,2) = gadd(p3,p4);
1680 : #else
1681 259979357 : if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1682 : { /* 3M formula */
1683 559356 : p3 = addii(xr,xi);
1684 559356 : p4 = addii(yr,yi);
1685 559356 : p1 = mulii(xr,yr);
1686 559356 : p2 = mulii(xi,yi);
1687 559356 : p3 = mulii(p3,p4);
1688 559356 : p4 = addii(p2,p1);
1689 559356 : tetpil = avma;
1690 559356 : gel(z,1) = subii(p1,p2);
1691 559356 : gel(z,2) = subii(p3,p4);
1692 559356 : if (!signe(gel(z,2)))
1693 113260 : return gc_INT((pari_sp)(z+3), gel(z,1));
1694 : }
1695 : else
1696 : { /* naive 4M formula: avoid all loss of accuracy */
1697 259420001 : p1 = gmul(xr,yr);
1698 259358970 : p2 = gmul(xi,yi);
1699 259352340 : p3 = gmul(xr,yi);
1700 259354875 : p4 = gmul(xi,yr);
1701 259363466 : tetpil = avma;
1702 259363466 : gel(z,1) = gsub(p1,p2);
1703 259225576 : gel(z,2) = gadd(p3,p4);
1704 259233932 : if (isintzero(gel(z,2)))
1705 : {
1706 50575 : cgiv(gel(z,2));
1707 50575 : return gc_upto((pari_sp)(z+3), gel(z,1));
1708 : }
1709 : }
1710 : #endif
1711 :
1712 259621629 : gc_slice_unsafe(av,tetpil, z+1,2); return z;
1713 : }
1714 : /* x,y PADIC */
1715 : static GEN
1716 17503064 : mulpp(GEN x, GEN y) {
1717 17503064 : GEN pd, p = padic_p(x), ux = padic_u(x), uy = padic_u(y);
1718 17503064 : long e = valp(x) + valp(y), dx, dy;
1719 :
1720 17503064 : if (!equalii(p, padic_p(y))) pari_err_OP("*",x,y);
1721 17502997 : if (!signe(ux) || !signe(uy)) return zeropadic(p, e);
1722 16930311 : dx = precp(x); dy = precp(y);
1723 16930311 : if (dx > dy) { pd = padic_pd(y); dx = dy; } else pd = padic_pd(x);
1724 16930311 : retmkpadic(Fp_mul(ux, uy, pd), icopy(p), icopy(pd), e, dx);
1725 : }
1726 : /* x,y QUAD, w^2 = - b w - c, where b = 0 or -1
1727 : * (ux + vx w)(uy + vy w) = ux uy - c vx vy + (ux vy + uy vx - b vx vy) w */
1728 : static GEN
1729 1127 : mulqq(GEN x, GEN y)
1730 : {
1731 1127 : GEN T = gel(x,1), b = gel(T,3), c = gel(T,2);
1732 1127 : GEN ux = gel(x,2), vx = gel(x,3), uy = gel(y,2), vy = gel(y,3);
1733 : GEN z, s, U, V, E, F;
1734 : pari_sp av, av2;
1735 1127 : z = cgetg(4, t_QUAD), gel(z,1) = ZX_copy(T); av = avma;
1736 1127 : U = gmul(ux, uy);
1737 1127 : V = gmul(vx, vy); s = gmul(c, V);
1738 1127 : E = gmul(gadd(ux, vx), gadd(uy, vy));
1739 1127 : F = signe(b)? U: gadd(U, V);
1740 : /* E - F = ux vy + uy vx - b vx vy */
1741 1127 : av2 = avma;
1742 1127 : gel(z,2) = gsub(U, s);
1743 1127 : gel(z,3) = gsub(E, F);
1744 1127 : gc_slice_unsafe(av, av2, z+2, 2); return z;
1745 : }
1746 :
1747 : GEN
1748 15655564 : mulcxI(GEN x)
1749 : {
1750 : GEN z;
1751 15655564 : switch(typ(x))
1752 : {
1753 1006784 : case t_INT:
1754 1006784 : if (!signe(x)) return gen_0;
1755 : case t_REAL: case t_FRAC:
1756 2633609 : return mkcomplex(gen_0, x);
1757 13021724 : case t_COMPLEX:
1758 13021724 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1759 13008275 : z = cgetg(3,t_COMPLEX);
1760 13008686 : gel(z,1) = gneg(gel(x,2));
1761 13009160 : gel(z,2) = gel(x,1); return z;
1762 56 : default:
1763 56 : return gmul(gen_I(), x);
1764 : }
1765 : }
1766 : GEN
1767 3079079 : mulcxmI(GEN x)
1768 : {
1769 : GEN z;
1770 3079079 : switch(typ(x))
1771 : {
1772 1092 : case t_INT:
1773 1092 : if (!signe(x)) return gen_0;
1774 : case t_REAL: case t_FRAC:
1775 117906 : return mkcomplex(gen_0, gneg(x));
1776 2758185 : case t_COMPLEX:
1777 2758185 : if (isintzero(gel(x,1))) return gel(x,2);
1778 2720939 : z = cgetg(3,t_COMPLEX);
1779 2720910 : gel(z,1) = gel(x,2);
1780 2720910 : gel(z,2) = gneg(gel(x,1)); return z;
1781 202988 : default:
1782 202988 : return gmul(mkcomplex(gen_0, gen_m1), x);
1783 : }
1784 : }
1785 : /* x * I^k */
1786 : GEN
1787 5527343 : mulcxpowIs(GEN x, long k)
1788 : {
1789 5527343 : switch (k & 3)
1790 : {
1791 1369256 : case 1: return mulcxI(x);
1792 1362290 : case 2: return gneg(x);
1793 1365891 : case 3: return mulcxmI(x);
1794 : }
1795 1429906 : return x;
1796 : }
1797 :
1798 : /* caller will assume l > 2 later */
1799 : static GEN
1800 5706368 : init_ser(long l, long v, long e)
1801 : {
1802 5706368 : GEN z = cgetg(l, t_SER);
1803 5706368 : z[1] = evalvalser(e) | evalvarn(v) | evalsigne(1); return z;
1804 : }
1805 :
1806 : /* fill in coefficients of t_SER z from coeffs of t_POL y */
1807 : static GEN
1808 5693993 : fill_ser(GEN z, GEN y)
1809 : {
1810 5693993 : long i, lx = lg(z), ly = lg(y); /* lx > 2 */
1811 5693993 : if (ly >= lx) {
1812 25570893 : for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
1813 : } else {
1814 335415 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1815 239573 : for ( ; i < lx; i++) gel(z,i) = gen_0;
1816 : }
1817 : /* dangerous case (already normalized), don't use normalizeser */
1818 5693993 : if (ly == 3 && !signe(y)) { setsigne(z, 0); return z; }
1819 5693895 : return normalizeser(z);
1820 : }
1821 : /* assume typ(x) = t_VEC */
1822 : static int
1823 98 : is_ext_qfr(GEN x)
1824 91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
1825 189 : && typ(gel(x,2)) == t_REAL; }
1826 :
1827 : GEN
1828 8823395922 : gmul(GEN x, GEN y)
1829 : {
1830 : long tx, ty, lx, ly, vx, vy, i, l;
1831 : pari_sp av, tetpil;
1832 : GEN z, p1;
1833 :
1834 8823395922 : if (x == y) return gsqr(x);
1835 7884690148 : tx = typ(x); ty = typ(y);
1836 7884690148 : if (tx == ty) switch(tx)
1837 : {
1838 3676353607 : case t_INT: return mulii(x,y);
1839 1893102882 : case t_REAL: return mulrr(x,y);
1840 1286534 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1841 1286534 : z = cgetg(3,t_INTMOD);
1842 1286534 : if (X==Y || equalii(X,Y))
1843 1286499 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1844 35 : gel(z,1) = gcdii(X,Y);
1845 35 : warn_coercion(X,Y,gel(z,1));
1846 35 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1847 35 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
1848 : }
1849 22792123 : case t_FRAC:
1850 : {
1851 22792123 : GEN x1 = gel(x,1), x2 = gel(x,2);
1852 22792123 : GEN y1 = gel(y,1), y2 = gel(y,2);
1853 22792123 : z=cgetg(3,t_FRAC);
1854 22792123 : p1 = gcdii(x1, y2);
1855 22792121 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1856 22792119 : p1 = gcdii(x2, y1);
1857 22792120 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1858 22792120 : tetpil = avma;
1859 22792120 : gel(z,2) = mulii(x2,y2);
1860 22792122 : gel(z,1) = mulii(x1,y1);
1861 22792121 : fix_frac_if_int_GC(z,tetpil); return z;
1862 : }
1863 269506822 : case t_COMPLEX: return mulcc(x, y);
1864 12509654 : case t_PADIC: return mulpp(x, y);
1865 882 : case t_QUAD:
1866 882 : if (!ZX_equal(gel(x,1), gel(y,1))) pari_err_OP("*",x,y);
1867 882 : return mulqq(x, y);
1868 13689031 : case t_FFELT: return FF_mul(x, y);
1869 10950787 : case t_POLMOD:
1870 10950787 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1871 8275842 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1872 2674945 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1873 134013680 : case t_POL:
1874 134013680 : vx = varn(x);
1875 134013680 : vy = varn(y);
1876 134013680 : if (vx != vy) {
1877 46850437 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1878 27724525 : else return RgX_Rg_mul(y, x);
1879 : }
1880 87163243 : return RgX_mul(x, y);
1881 :
1882 4411917 : case t_SER: {
1883 4411917 : vx = varn(x);
1884 4411917 : vy = varn(y);
1885 4411917 : if (vx != vy) {
1886 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1887 3675 : return mul_ser_scal(y, x);
1888 : }
1889 4408242 : lx = minss(lg(x), lg(y));
1890 4408242 : if (lx == 2) return zeroser(vx, valser(x)+valser(y));
1891 4408228 : av = avma; z = init_ser(lx, vx, valser(x)+valser(y));
1892 4408228 : x = ser2pol_i(x, lx);
1893 4408228 : y = ser2pol_i(y, lx);
1894 4408228 : y = RgXn_mul(x, y, lx-2);
1895 4408228 : return gc_GEN(av, fill_ser(z,y));
1896 : }
1897 42 : case t_VEC:
1898 42 : if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
1899 : /* fall through, handle extended t_QFB */
1900 198 : case t_QFB: return qfbcomp(x,y);
1901 6723048 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1902 4125769 : case t_MAT: return RgM_mul(x, y);
1903 :
1904 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1905 1631 : z = cgetg_copy(x, &l);
1906 1631 : if (l != lg(y)) break;
1907 17325 : for (i=1; i<l; i++)
1908 : {
1909 15694 : long yi = y[i];
1910 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1911 15694 : z[i] = x[yi];
1912 : }
1913 1631 : return z;
1914 :
1915 0 : default:
1916 0 : pari_err_TYPE2("*",x,y);
1917 : }
1918 : /* tx != ty */
1919 1837272264 : if (is_const_t(ty) && is_const_t(tx)) {
1920 1670981498 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1921 1670981498 : switch(tx) {
1922 1541494334 : case t_INT:
1923 : switch(ty)
1924 : {
1925 994863894 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1926 1342479 : case t_INTMOD:
1927 1342479 : z = cgetg(3, t_INTMOD);
1928 1342479 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1929 106447736 : case t_FRAC:
1930 106447736 : if (!signe(x)) return gen_0;
1931 89312291 : z=cgetg(3,t_FRAC);
1932 89312471 : p1 = gcdii(x,gel(y,2));
1933 89311627 : if (equali1(p1))
1934 : {
1935 49455130 : set_avma((pari_sp)z);
1936 49455113 : gel(z,2) = icopy(gel(y,2));
1937 49455188 : gel(z,1) = mulii(gel(y,1), x);
1938 : }
1939 : else
1940 : {
1941 39856794 : x = diviiexact(x,p1); tetpil = avma;
1942 39856196 : gel(z,2) = diviiexact(gel(y,2), p1);
1943 39855927 : gel(z,1) = mulii(gel(y,1), x);
1944 39856482 : fix_frac_if_int_GC(z,tetpil);
1945 : }
1946 89312226 : return z;
1947 436070926 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1948 4198901 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1949 1904 : case t_QUAD: return mulRq(x,y);
1950 1607079 : case t_FFELT: return FF_Z_mul(y,x);
1951 : }
1952 :
1953 : case t_REAL:
1954 121704974 : switch(ty)
1955 : {
1956 37790599 : case t_FRAC: return mulrfrac(x, y);
1957 83914333 : case t_COMPLEX: return mulRc(x, y);
1958 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1959 21 : default: pari_err_TYPE2("*",x,y);
1960 : }
1961 :
1962 8036 : case t_INTMOD:
1963 : switch(ty)
1964 : {
1965 7133 : case t_FRAC: { GEN X = gel(x,1);
1966 7133 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1967 7133 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1968 : }
1969 49 : case t_COMPLEX: return mulRc(x,y);
1970 259 : case t_PADIC: { GEN X = gel(x,1);
1971 259 : z = cgetg(3, t_INTMOD);
1972 259 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1973 : }
1974 63 : case t_QUAD: return mulRq(x, y);
1975 532 : case t_FFELT:
1976 532 : if (!equalii(gel(x,1),FF_p_i(y)))
1977 0 : pari_err_OP("*",x,y);
1978 532 : return FF_Z_mul(y,gel(x,2));
1979 : }
1980 :
1981 : case t_FRAC:
1982 : switch(ty)
1983 : {
1984 5375788 : case t_COMPLEX: return mulRc(x, y);
1985 2582752 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
1986 343 : case t_QUAD: return mulRq(x, y);
1987 2079 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
1988 : }
1989 :
1990 : case t_FFELT:
1991 35 : pari_err_TYPE2("*",x,y);
1992 :
1993 21 : case t_COMPLEX:
1994 : switch(ty)
1995 : {
1996 14 : case t_PADIC:
1997 14 : return Zp_nosquare_m1(padic_p(y))? mulRc(y, x): mulTp(x, y);
1998 7 : case t_QUAD:
1999 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
2000 7 : return mulqf(y, x, lx);
2001 : }
2002 :
2003 : case t_PADIC: /* ty == t_QUAD */
2004 28 : return (kro_quad(y,padic_p(x))== -1)? mulRq(x, y): mulTp(y, x);
2005 : }
2006 : }
2007 :
2008 167981095 : if (is_matvec_t(ty)) switch(tx)
2009 : {
2010 122011 : case t_VEC:
2011 : switch(ty) {
2012 14826 : case t_COL: return RgV_RgC_mul(x,y);
2013 107185 : case t_MAT: return RgV_RgM_mul(x,y);
2014 : }
2015 0 : break;
2016 1687 : case t_COL:
2017 : switch(ty) {
2018 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2019 0 : case t_MAT: return RgC_RgM_mul(x,y);
2020 : }
2021 0 : break;
2022 880392 : case t_MAT:
2023 : switch(ty) {
2024 0 : case t_VEC: return RgM_RgV_mul(x,y);
2025 880392 : case t_COL: return RgM_RgC_mul(x,y);
2026 : }
2027 : default:
2028 5419238 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
2029 5419244 : return mul_self_scal(y, x);
2030 : }
2031 165020531 : if (is_matvec_t(tx))
2032 : {
2033 2126413 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2034 2126496 : return mul_self_scal(x, y);
2035 : }
2036 162894150 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2037 : /* tx < ty, !ismatvec(x and y) */
2038 :
2039 162894150 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2040 3215290 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2041 159678860 : if (is_scalar_t(tx)) {
2042 156209705 : if (tx == t_POLMOD) {
2043 3159530 : vx = varn(gel(x,1));
2044 3159530 : vy = gvar(y);
2045 3159530 : if (vx != vy) {
2046 2907313 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2047 23093 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2048 : }
2049 : /* error if ty == t_SER */
2050 252217 : av = avma; y = gmod(y, gel(x,1));
2051 252210 : return gc_upto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2052 : }
2053 153050175 : return mul_scal(y, x, ty);
2054 : }
2055 :
2056 : /* x and y are not scalars, nor matvec */
2057 3469128 : vx = gvar(x);
2058 3469155 : vy = gvar(y);
2059 3469155 : if (vx != vy) /* x or y is treated as a scalar */
2060 2785619 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2061 2785619 : : mul_scal(y, x, ty);
2062 : /* vx = vy */
2063 2045731 : switch(tx)
2064 : {
2065 2045696 : case t_POL:
2066 : switch (ty)
2067 : {
2068 6741 : case t_SER:
2069 : {
2070 : long v;
2071 6741 : av = avma; v = RgX_valrem(x, &x);
2072 6741 : if (v == LONG_MAX) return gc_upto(av, Rg_get_0(x));
2073 6727 : v += valser(y); ly = lg(y);
2074 6727 : if (ly == 2) { set_avma(av); return zeroser(vx, v); }
2075 6496 : if (degpol(x))
2076 : {
2077 2030 : z = init_ser(ly, vx, v);
2078 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2079 2030 : return gc_GEN(av, fill_ser(z, x));
2080 : }
2081 : /* take advantage of x = c*t^v */
2082 4466 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2083 4466 : setvalser(y, v); return y;
2084 : }
2085 :
2086 2038955 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2087 : }
2088 0 : break;
2089 :
2090 35 : case t_SER:
2091 : switch (ty)
2092 : {
2093 35 : case t_RFRAC:
2094 35 : av = avma;
2095 35 : return gc_upto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2096 : }
2097 0 : break;
2098 : }
2099 0 : pari_err_TYPE2("*",x,y);
2100 : return NULL; /* LCOV_EXCL_LINE */
2101 : }
2102 :
2103 : /* return a nonnormalized result */
2104 : GEN
2105 112759 : sqr_ser_part(GEN x, long l1, long l2)
2106 : {
2107 : long i, j, l;
2108 : pari_sp av;
2109 : GEN Z, z, p1, p2;
2110 : long mi;
2111 112759 : if (l2 < l1) return zeroser(varn(x), 2*valser(x));
2112 112745 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2113 112745 : Z = cgetg(l2-l1+3, t_SER);
2114 112745 : Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
2115 112745 : z = Z + 2-l1;
2116 112745 : x += 2; mi = 0;
2117 425418 : for (i=0; i<l1; i++)
2118 : {
2119 312673 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2120 : }
2121 :
2122 720520 : for (i=l1; i<=l2; i++)
2123 : {
2124 607775 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2125 607775 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2126 1256988 : for (j=i-mi; j<=minss(l,mi); j++)
2127 649213 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2128 607775 : p1 = gshift(p1,1);
2129 607775 : if ((i&1) == 0 && p2[i>>1])
2130 93669 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2131 607775 : gel(z,i) = gc_upto(av,p1);
2132 : }
2133 112745 : return Z;
2134 : }
2135 :
2136 : /* (u + v X)^2 mod (X^2 + bX + c), b = 0 or -1
2137 : * = u^2 - c v^2 + (2uv - b v^2) X */
2138 : static GEN
2139 70 : sqrq(GEN x)
2140 : {
2141 70 : GEN T = gel(x,1), c = gel(T,2), b = gel(T,3);
2142 70 : GEN u = gel(x,2), v = gel(x,3), U, V, E, F, s, z;
2143 : pari_sp av, av2;
2144 :
2145 70 : z = cgetg(4, t_QUAD), gel(z,1) = ZX_copy(T); av = avma;
2146 70 : U = gsqr(u);
2147 70 : V = gsqr(v); s = gmul(c, V);
2148 70 : E = gmul(u, v);
2149 70 : F = signe(b)? gadd(E, V): E; /* u v - b v^2 */
2150 70 : av2 = avma;
2151 70 : gel(z,2) = gsub(U, s);
2152 70 : gel(z,3) = gadd(E, F);
2153 70 : gc_slice_unsafe(av, av2, z+2, 2); return z;
2154 : }
2155 :
2156 : GEN
2157 1181236865 : gsqr(GEN x)
2158 : {
2159 : long i, lx;
2160 : pari_sp av, tetpil;
2161 : GEN z, p1, p2, p3;
2162 :
2163 1181236865 : switch(typ(x))
2164 : {
2165 968089443 : case t_INT: return sqri(x);
2166 195827170 : case t_REAL: return sqrr(x);
2167 142467 : case t_INTMOD: { GEN X = gel(x,1);
2168 142467 : z = cgetg(3,t_INTMOD);
2169 142467 : gel(z,2) = gc_INT((pari_sp)z, remii(sqri(gel(x,2)), X));
2170 142465 : gel(z,1) = icopy(X); return z;
2171 : }
2172 4143033 : case t_FRAC: return sqrfrac(x);
2173 :
2174 8282013 : case t_COMPLEX:
2175 8282013 : if (isintzero(gel(x,1))) {
2176 243286 : av = avma;
2177 243286 : return gc_upto(av, gneg(gsqr(gel(x,2))));
2178 : }
2179 8038719 : z = cgetg(3,t_COMPLEX); av = avma;
2180 8038697 : p1 = gadd(gel(x,1),gel(x,2));
2181 8038462 : p2 = gsub(gel(x,1), gel(x,2));
2182 8038491 : p3 = gmul(gel(x,1),gel(x,2));
2183 8038612 : tetpil = avma;
2184 8038612 : gel(z,1) = gmul(p1,p2);
2185 8038692 : gel(z,2) = gshift(p3,1); gc_slice_unsafe(av,tetpil,z+1,2); return z;
2186 :
2187 4795 : case t_PADIC:
2188 : {
2189 4795 : GEN u = padic_u(x), p = padic_p(x), pd = padic_pd(x);
2190 4795 : long v2 = 2*valp(x), d = precp(x);
2191 4795 : if (!signe(u)) { x = gcopy(x); setvalp(x, v2); return x; }
2192 4739 : if (!absequaliu(p, 2))
2193 3206 : retmkpadic(Fp_sqr(u, pd), icopy(p), icopy(pd), v2, d);
2194 : /* p = 2*/
2195 1533 : if (d == 1) /* (1 + O(2))^2 = 1 + O(2^3) */
2196 7 : retmkpadic(gen_1, gen_2, utoipos(8), v2, 3);
2197 1526 : retmkpadic(remi2n(sqri(u), d + 1), gen_2, int2n(d + 1), v2, d + 1);
2198 : }
2199 70 : case t_QUAD: return sqrq(x);
2200 :
2201 419385 : case t_POLMOD: return sqr_polmod(gel(x,1), gel(x,2));
2202 :
2203 2967119 : case t_FFELT: return FF_sqr(x);
2204 :
2205 1351465 : case t_POL: return RgX_sqr(x);
2206 :
2207 35042 : case t_SER:
2208 35042 : lx = lg(x);
2209 35042 : if (ser_isexactzero(x)) {
2210 21 : GEN z = gcopy(x);
2211 21 : setvalser(z, 2*valser(x));
2212 21 : return z;
2213 : }
2214 35021 : if (lx < 40)
2215 34762 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2216 : else
2217 : {
2218 259 : pari_sp av = avma;
2219 259 : GEN z = init_ser(lx, varn(x), 2*valser(x));
2220 259 : x = ser2pol_i(x, lx);
2221 259 : x = RgXn_sqr(x, lx-2);
2222 259 : return gc_GEN(av, fill_ser(z,x));
2223 : }
2224 :
2225 63 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2226 63 : gel(z,1) = gsqr(gel(x,1));
2227 63 : gel(z,2) = gsqr(gel(x,2)); return z;
2228 :
2229 1169 : case t_MAT: return RgM_sqr(x);
2230 28 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
2231 : /* fall through handle extended t_QFB */
2232 36 : case t_QFB: return qfbsqr(x);
2233 658 : case t_VECSMALL:
2234 658 : z = cgetg_copy(x, &lx);
2235 16289 : for (i=1; i<lx; i++)
2236 : {
2237 15631 : long xi = x[i];
2238 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2239 15631 : z[i] = x[xi];
2240 : }
2241 658 : return z;
2242 : }
2243 7 : pari_err_TYPE2("*",x,x);
2244 : return NULL; /* LCOV_EXCL_LINE */
2245 : }
2246 :
2247 : /********************************************************************/
2248 : /** **/
2249 : /** DIVISION **/
2250 : /** **/
2251 : /********************************************************************/
2252 : static GEN
2253 32139 : div_rfrac_scal(GEN x, GEN y)
2254 : {
2255 32139 : pari_sp av = avma;
2256 32139 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2257 32139 : return gc_upto(av, gred_rfrac_simple(gel(x,1), d));
2258 : }
2259 : static GEN
2260 37767 : div_scal_rfrac(GEN x, GEN y)
2261 : {
2262 37767 : GEN y1 = gel(y,1), y2 = gel(y,2);
2263 37767 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2264 : {
2265 14 : if (degpol(y1))
2266 : {
2267 14 : pari_sp av = avma;
2268 14 : GEN _1 = Rg_get_1(x);
2269 14 : if (transtype(_1)) y1 = gmul(y1, _1);
2270 14 : return gc_upto(av, gred_rfrac_simple(gmul(x, y2), y1));
2271 : }
2272 0 : y1 = gel(y1,2);
2273 : }
2274 37753 : return RgX_Rg_mul(y2, gdiv(x,y1));
2275 : }
2276 : static GEN
2277 1187237 : div_rfrac(GEN x, GEN y)
2278 1187237 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2279 :
2280 : /* x != 0 */
2281 : static GEN
2282 1338593 : div_ser_scal(GEN x, GEN t)
2283 : {
2284 1338593 : if (ser_isexactzero(x))
2285 : {
2286 28 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gdiv(gel(x,2), t), varn(x), 1);
2287 28 : setvalser(z, valser(x)); return z;
2288 : }
2289 6202113 : pari_APPLY_ser(gdiv(gel(x,i), t));
2290 : }
2291 : GEN
2292 7658 : ser_normalize(GEN x)
2293 : {
2294 7658 : long i, lx = lg(x);
2295 : GEN c, z;
2296 7658 : if (lx == 2) return x;
2297 7658 : c = gel(x,2); if (gequal1(c)) return x;
2298 7581 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2299 108752 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2300 7581 : return z;
2301 : }
2302 :
2303 : /* y != 0 */
2304 : static GEN
2305 7122689 : div_T_scal(GEN x, GEN y, long tx) {
2306 7122689 : switch(tx)
2307 : {
2308 5755807 : case t_POL: return RgX_Rg_div(x, y);
2309 1338586 : case t_SER: return div_ser_scal(x, y);
2310 28296 : case t_RFRAC: return div_rfrac_scal(x,y);
2311 : }
2312 0 : pari_err_TYPE2("/",x,y);
2313 : return NULL; /* LCOV_EXCL_LINE */
2314 : }
2315 :
2316 : static GEN
2317 9257030 : div_scal_pol(GEN x, GEN y) {
2318 9257030 : long ly = lg(y);
2319 : pari_sp av;
2320 : GEN _1;
2321 9257030 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2322 9178160 : if (isrationalzero(x)) return zeropol(varn(y));
2323 7131376 : av = avma;
2324 7131376 : _1 = Rg_get_1(x); if (transtype(_1)) y = gmul(y, _1);
2325 7131376 : return gc_upto(av, gred_rfrac_simple(x,y));
2326 : }
2327 : static GEN
2328 18550 : div_scal_ser(GEN x, GEN y)
2329 : {
2330 18550 : pari_sp av = avma;
2331 18550 : GEN _1 = Rg_get_1(x);
2332 18550 : if (transtype(_1)) y = gmul(y, _1);
2333 18550 : return gc_upto(av, gmul(x, ser_inv(y)));
2334 : }
2335 : static GEN
2336 9264727 : div_scal_T(GEN x, GEN y, long ty) {
2337 9264727 : switch(ty)
2338 : {
2339 9209173 : case t_POL: return div_scal_pol(x, y);
2340 18543 : case t_SER: return div_scal_ser(x, y);
2341 37011 : case t_RFRAC: return div_scal_rfrac(x, y);
2342 : }
2343 0 : pari_err_TYPE2("/",x,y);
2344 : return NULL; /* LCOV_EXCL_LINE */
2345 : }
2346 :
2347 : /* assume tx = ty = t_SER, same variable vx */
2348 : static GEN
2349 771538 : div_ser(GEN x, GEN y, long vx)
2350 : {
2351 771538 : long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
2352 771538 : GEN y0 = y, z;
2353 771538 : pari_sp av = avma;
2354 :
2355 771538 : if (!signe(y)) pari_err_INV("div_ser", y);
2356 771531 : if (ser_isexactzero(x))
2357 : {
2358 59892 : if (lx == 2) return zeroser(vx, v);
2359 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2360 147 : setvalser(z, v); return z;
2361 : }
2362 711639 : if (lx < ly) ly = lx;
2363 711639 : y = ser2pol_i_normalize(y, ly, &e);
2364 711639 : if (e)
2365 : {
2366 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2367 7 : ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
2368 : }
2369 711639 : z = init_ser(ly, vx, v);
2370 711639 : if (ly == 3)
2371 : {
2372 12375 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2373 12375 : if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
2374 12375 : return gc_upto(av, z);
2375 : }
2376 699264 : x = ser2pol_i(x, ly);
2377 699264 : y = RgXn_div_i(x, y, ly-2);
2378 699264 : return gc_GEN(av, fill_ser(z,y));
2379 : }
2380 : /* x,y compatible PADIC */
2381 : static GEN
2382 13212381 : divpp(GEN x, GEN y)
2383 : {
2384 13212381 : GEN M, ux = padic_u(x), uy = padic_u(y), p = padic_p(x);
2385 13212381 : long a = precp(x), b = precp(y), v = valp(x) - valp(y);
2386 :
2387 13212381 : if (!signe(uy)) pari_err_INV("divpp",y);
2388 13212397 : if (!signe(ux)) return zeropadic(p, v);
2389 13212068 : if (a > b) M = padic_pd(y); else { M = padic_pd(x); b = a; }
2390 13212068 : retmkpadic(Fp_div(ux, uy, M), icopy(padic_p(x)), icopy(M), v, b);
2391 : }
2392 : static GEN
2393 36764 : div_polmod_same(GEN T, GEN x, GEN y)
2394 : {
2395 36764 : long v = varn(T);
2396 36764 : GEN a, z = cgetg(3, t_POLMOD);
2397 36764 : gel(z,1) = RgX_copy(T);
2398 36764 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2399 11102 : a = gdiv(x, y);
2400 25662 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2401 105 : {
2402 105 : pari_sp av = avma;
2403 105 : a = gc_upto(av, gmul(x, RgXQ_inv(y, T)));
2404 : }
2405 25557 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2406 5278 : {
2407 5278 : pari_sp av = avma;
2408 5278 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2409 5278 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2410 5278 : a = gc_upto(av, a);
2411 : }
2412 : else
2413 : {
2414 20279 : pari_sp av = avma;
2415 20279 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2416 20279 : a = gc_upto(av, a);
2417 : }
2418 36764 : gel(z,2) = a; return z;
2419 : }
2420 : GEN
2421 438379270 : gdiv(GEN x, GEN y)
2422 : {
2423 438379270 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2424 : pari_sp av, tetpil;
2425 : GEN z, p1;
2426 :
2427 438379270 : if (tx == ty) switch(tx)
2428 : {
2429 91170641 : case t_INT:
2430 91170641 : return Qdivii(x,y);
2431 :
2432 94388079 : case t_REAL: return divrr(x,y);
2433 25719 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2434 25719 : z = cgetg(3,t_INTMOD);
2435 25719 : if (X==Y || equalii(X,Y))
2436 25705 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2437 14 : gel(z,1) = gcdii(X,Y);
2438 14 : warn_coercion(X,Y,gel(z,1));
2439 14 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2440 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
2441 : }
2442 4042526 : case t_FRAC: {
2443 4042526 : GEN x1 = gel(x,1), x2 = gel(x,2);
2444 4042526 : GEN y1 = gel(y,1), y2 = gel(y,2);
2445 4042526 : z = cgetg(3, t_FRAC);
2446 4042526 : p1 = gcdii(x1, y1);
2447 4042524 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2448 4042523 : p1 = gcdii(x2, y2);
2449 4042523 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2450 4042519 : tetpil = avma;
2451 4042519 : gel(z,2) = mulii(x2,y1);
2452 4042521 : gel(z,1) = mulii(x1,y2);
2453 4042523 : normalize_frac(z);
2454 4042523 : fix_frac_if_int_GC(z,tetpil);
2455 4042526 : return z;
2456 : }
2457 9479745 : case t_COMPLEX:
2458 9479745 : if (isintzero(gel(y,1)))
2459 : {
2460 197113 : y = gel(y,2);
2461 197113 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2462 116039 : z = cgetg(3,t_COMPLEX);
2463 116039 : gel(z,1) = gdiv(gel(x,2), y);
2464 116039 : av = avma;
2465 116039 : gel(z,2) = gc_upto(av, gneg(gdiv(gel(x,1), y)));
2466 116039 : return z;
2467 : }
2468 9282629 : av = avma;
2469 9282629 : return gc_upto(av, gdiv(mulcc(x, conj_i(y)), cxnorm(y)));
2470 :
2471 587462 : case t_PADIC:
2472 587462 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("/",x,y);
2473 587462 : return divpp(x, y);
2474 :
2475 259 : case t_QUAD:
2476 259 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2477 245 : av = avma;
2478 245 : return gc_upto(av, gdiv(mulqq(x, conj_i(y)), quadnorm(y)));
2479 :
2480 246054 : case t_FFELT: return FF_div(x,y);
2481 :
2482 41139 : case t_POLMOD:
2483 41139 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2484 36764 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2485 : else {
2486 4375 : av = avma;
2487 4375 : z = gc_upto(av, gmul(x, ginv(y)));
2488 : }
2489 41139 : return z;
2490 :
2491 18571576 : case t_POL:
2492 18571576 : vx = varn(x);
2493 18571576 : vy = varn(y);
2494 18571576 : if (vx != vy) {
2495 154399 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2496 47857 : else return div_scal_pol(x, y);
2497 : }
2498 18417177 : if (!signe(y)) pari_err_INV("gdiv",y);
2499 18417177 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2500 18288207 : av = avma;
2501 18288207 : return gc_upto(av, gred_rfrac2(x,y));
2502 :
2503 771559 : case t_SER:
2504 771559 : vx = varn(x);
2505 771559 : vy = varn(y);
2506 771559 : if (vx != vy) {
2507 21 : if (varncmp(vx, vy) < 0)
2508 : {
2509 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2510 7 : return div_ser_scal(x, y);
2511 : }
2512 7 : return div_scal_ser(x, y);
2513 : }
2514 771538 : return div_ser(x, y, vx);
2515 1191696 : case t_RFRAC:
2516 1191696 : vx = varn(gel(x,2));
2517 1191696 : vy = varn(gel(y,2));
2518 1191696 : if (vx != vy) {
2519 4459 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2520 756 : else return div_scal_rfrac(x, y);
2521 : }
2522 1187237 : return div_rfrac(x,y);
2523 :
2524 21 : case t_VEC: /* handle extended t_QFB */
2525 21 : case t_QFB: av = avma; return gc_upto(av, qfbcomp(x, ginv(y)));
2526 :
2527 28 : case t_MAT:
2528 28 : p1 = RgM_div(x,y);
2529 28 : if (!p1) pari_err_INV("gdiv",y);
2530 21 : return p1;
2531 :
2532 0 : default: pari_err_TYPE2("/",x,y);
2533 : }
2534 :
2535 217867869 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2536 : {
2537 3801609 : long s = signe(x);
2538 3801609 : if (!s) {
2539 517887 : if (gequal0(y)) pari_err_INV("gdiv",y);
2540 517887 : switch (ty)
2541 : {
2542 514702 : default: return gen_0;
2543 70 : case t_INTMOD:
2544 70 : z = cgetg(3,t_INTMOD);
2545 70 : gel(z,1) = icopy(gel(y,1));
2546 70 : gel(z,2) = gen_0; return z;
2547 3115 : case t_FFELT: return FF_zero(y);
2548 : }
2549 : }
2550 3283722 : if (is_pm1(x)) {
2551 1325163 : if (s > 0) return ginv(y);
2552 228412 : av = avma; return gc_upto(av, ginv(gneg(y)));
2553 : }
2554 1958559 : switch(ty)
2555 : {
2556 644410 : case t_REAL: return divir(x,y);
2557 21 : case t_INTMOD:
2558 21 : z = cgetg(3, t_INTMOD);
2559 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2560 792284 : case t_FRAC:
2561 792284 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2562 792284 : if (equali1(p1))
2563 : {
2564 271463 : set_avma((pari_sp)z);
2565 271463 : gel(z,2) = icopy(gel(y,1));
2566 271463 : gel(z,1) = mulii(gel(y,2), x);
2567 271463 : normalize_frac(z);
2568 271463 : fix_frac_if_int(z);
2569 : }
2570 : else
2571 : {
2572 520821 : x = diviiexact(x,p1); tetpil = avma;
2573 520821 : gel(z,2) = diviiexact(gel(y,1), p1);
2574 520821 : gel(z,1) = mulii(gel(y,2), x);
2575 520821 : normalize_frac(z);
2576 520821 : fix_frac_if_int_GC(z,tetpil);
2577 : }
2578 792284 : return z;
2579 :
2580 420 : case t_FFELT: return Z_FF_div(x,y);
2581 519693 : case t_COMPLEX: return divRc(x,y);
2582 1505 : case t_PADIC: return divTp(x, y);
2583 231 : case t_QUAD: return divRq(x,y);
2584 : }
2585 : }
2586 214066254 : if (gequal0(y))
2587 : {
2588 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2589 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2590 : }
2591 :
2592 214075936 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2593 : {
2594 173421087 : case t_REAL:
2595 173421087 : switch(ty)
2596 : {
2597 171053458 : case t_INT: return divri(x,y);
2598 522692 : case t_FRAC:
2599 522692 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2600 522692 : return gc_leaf(av, z);
2601 1844984 : case t_COMPLEX: return divRc(x, y);
2602 42 : case t_QUAD: return divfq(x, y, realprec(x));
2603 18 : default: pari_err_TYPE2("/",x,y);
2604 : }
2605 :
2606 273 : case t_INTMOD:
2607 : switch(ty)
2608 : {
2609 189 : case t_INT:
2610 189 : z = cgetg(3, t_INTMOD);
2611 189 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2612 28 : case t_FRAC: { GEN X = gel(x,1);
2613 28 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2614 28 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2615 : }
2616 7 : case t_FFELT:
2617 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2618 0 : pari_err_OP("/",x,y);
2619 7 : return Z_FF_div(gel(x,2),y);
2620 :
2621 0 : case t_COMPLEX: return divRc(x,y);
2622 14 : case t_QUAD: return divRq(x,y);
2623 :
2624 7 : case t_PADIC: { GEN X = gel(x,1);
2625 7 : z = cgetg(3, t_INTMOD);
2626 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2627 : }
2628 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2629 : }
2630 :
2631 : case t_FRAC:
2632 : switch(ty)
2633 : {
2634 2102263 : case t_INT: z = cgetg(3, t_FRAC);
2635 2102263 : p1 = gcdii(y,gel(x,1));
2636 2102263 : if (equali1(p1))
2637 : {
2638 832574 : set_avma((pari_sp)z); tetpil = 0;
2639 832574 : gel(z,1) = icopy(gel(x,1));
2640 : }
2641 : else
2642 : {
2643 1269689 : y = diviiexact(y,p1); tetpil = avma;
2644 1269689 : gel(z,1) = diviiexact(gel(x,1), p1);
2645 : }
2646 2102263 : gel(z,2) = mulii(gel(x,2),y);
2647 2102263 : normalize_frac(z);
2648 2102263 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2649 2102263 : return z;
2650 :
2651 59572 : case t_REAL:
2652 59572 : av = avma;
2653 59572 : return gc_leaf(av, divir(gel(x,1), mulri(y,gel(x,2))));
2654 :
2655 7 : case t_INTMOD: { GEN Y = gel(y,1);
2656 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2657 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2658 : }
2659 :
2660 28 : case t_FFELT: av=avma;
2661 28 : return gc_upto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2662 :
2663 882 : case t_COMPLEX: return divRc(x, y);
2664 :
2665 2141 : case t_PADIC:
2666 2141 : if (!signe(gel(x,1))) return gen_0;
2667 2141 : return divTp(x, y);
2668 :
2669 42 : case t_QUAD: return divRq(x, y);
2670 : }
2671 :
2672 : case t_FFELT:
2673 161 : switch (ty)
2674 : {
2675 77 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2676 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2677 7 : case t_INTMOD:
2678 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2679 0 : pari_err_OP("/",x,y);
2680 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2681 49 : default:
2682 49 : pari_err_TYPE2("/",x,y);
2683 : }
2684 0 : break;
2685 :
2686 13780755 : case t_COMPLEX:
2687 : switch(ty)
2688 : {
2689 13780758 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2690 0 : case t_INTMOD: return mulRc(ginv(y), x);
2691 0 : case t_PADIC:
2692 0 : return Zp_nosquare_m1(padic_p(y))? divcR(x,y): divTp(x, y);
2693 0 : case t_QUAD:
2694 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2695 0 : return divfq(x, y, lx);
2696 : }
2697 :
2698 : case t_PADIC:
2699 : switch(ty)
2700 : {
2701 1205008 : case t_INT: case t_FRAC: { GEN p = padic_p(x);
2702 1204917 : return signe(padic_u(x))? divpT(x, y)
2703 2409925 : : zeropadic(p, valp(x) - Q_pval(y,p));
2704 : }
2705 7 : case t_INTMOD: { GEN Y = gel(y,1);
2706 7 : z = cgetg(3, t_INTMOD);
2707 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2708 : }
2709 0 : case t_COMPLEX: return divRc(x,y);
2710 14 : case t_QUAD: return divRq(x,y);
2711 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2712 : }
2713 :
2714 : case t_QUAD:
2715 : switch (ty)
2716 : {
2717 1204 : case t_INT: case t_INTMOD: case t_FRAC:
2718 1204 : z = cgetg(4,t_QUAD);
2719 1204 : gel(z,1) = ZX_copy(gel(x,1));
2720 1204 : gel(z,2) = gdiv(gel(x,2), y);
2721 1204 : gel(z,3) = gdiv(gel(x,3), y); return z;
2722 63 : case t_REAL: return divqf(x, y, realprec(y));
2723 28 : case t_PADIC: return divTp(x, y);
2724 0 : case t_COMPLEX:
2725 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2726 0 : return divqf(x, y, ly);
2727 : }
2728 : }
2729 23502228 : switch(ty) {
2730 582487 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2731 582487 : av = avma; return gc_upto(av, gmul(x, ginv(y)));
2732 42 : case t_MAT:
2733 42 : av = avma; p1 = RgM_inv(y);
2734 35 : if (!p1) pari_err_INV("gdiv",y);
2735 28 : return gc_upto(av, gmul(x, p1));
2736 0 : case t_VEC: case t_COL:
2737 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2738 0 : pari_err_TYPE2("/",x,y);
2739 : }
2740 22919699 : switch(tx) {
2741 4422582 : case t_VEC: case t_COL: case t_MAT:
2742 14369441 : pari_APPLY_same(gdiv(gel(x,i),y));
2743 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2744 0 : pari_err_TYPE2("/",x,y);
2745 : }
2746 :
2747 18497117 : vy = gvar(y);
2748 18498003 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2749 218133 : vx = varn(X);
2750 218133 : if (vx != vy) {
2751 217363 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2752 216950 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2753 : }
2754 : /* y is POL, SER or RFRAC */
2755 770 : av = avma;
2756 770 : switch(ty)
2757 : {
2758 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2759 770 : default: y = ginvmod(gmod(y,X), X);
2760 : }
2761 763 : return gc_upto(av, mul_polmod_same(X, gel(x,2), y));
2762 : }
2763 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2764 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2765 : * variable NO_VARIABLE, then the operation is incorrect */
2766 18279870 : vx = gvar(x);
2767 18279866 : if (vx != vy) { /* includes cases where one is scalar */
2768 16387006 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2769 9264313 : else return div_scal_T(x, y, ty);
2770 : }
2771 1892860 : switch(tx)
2772 : {
2773 1046171 : case t_POL:
2774 : switch(ty)
2775 : {
2776 28 : case t_SER:
2777 : {
2778 28 : GEN y0 = y;
2779 : long v;
2780 28 : av = avma; v = RgX_valrem(x, &x);
2781 28 : if (v == LONG_MAX) return gc_upto(av, Rg_get_0(x));
2782 14 : v -= valser(y); ly = lg(y); /* > 2 */
2783 14 : y = ser2pol_i_normalize(y, ly, &i);
2784 14 : if (i)
2785 : {
2786 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2787 7 : ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
2788 : }
2789 7 : z = init_ser(ly, vx, v);
2790 7 : return gc_GEN(av, fill_ser(z, RgXn_div(x, y, ly-2)));
2791 : }
2792 :
2793 1046143 : case t_RFRAC:
2794 : {
2795 1046143 : GEN y1 = gel(y,1), y2 = gel(y,2);
2796 1046143 : if (typ(y1) == t_POL && varn(y1) == vx)
2797 1143 : return mul_rfrac_scal(y2, y1, x);
2798 1045000 : av = avma;
2799 1045000 : return gc_upto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2800 : }
2801 : }
2802 0 : break;
2803 :
2804 584324 : case t_SER:
2805 : switch(ty)
2806 : {
2807 584310 : case t_POL:
2808 : {
2809 584310 : long v = valser(x);
2810 584310 : lx = lg(x);
2811 584310 : if (lx == 2) return zeroser(vx, v - RgX_val(y));
2812 584205 : av = avma;
2813 584205 : x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
2814 584205 : z = init_ser(lx, vx, v);
2815 584205 : if (!signe(x)) setsigne(z,0);
2816 584205 : return gc_GEN(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
2817 : }
2818 14 : case t_RFRAC:
2819 14 : av = avma;
2820 14 : return gc_upto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2821 : }
2822 0 : break;
2823 :
2824 262344 : case t_RFRAC:
2825 : switch(ty)
2826 : {
2827 262344 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2828 0 : case t_SER:
2829 0 : av = avma;
2830 0 : return gc_upto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2831 : }
2832 0 : break;
2833 : }
2834 21 : pari_err_TYPE2("/",x,y);
2835 : return NULL; /* LCOV_EXCL_LINE */
2836 : }
2837 :
2838 : /********************************************************************/
2839 : /** **/
2840 : /** SIMPLE MULTIPLICATION **/
2841 : /** **/
2842 : /********************************************************************/
2843 : GEN
2844 231771818 : gmulsg(long s, GEN x)
2845 : {
2846 : pari_sp av;
2847 : long i;
2848 : GEN z;
2849 :
2850 231771818 : switch(typ(x))
2851 : {
2852 137647425 : case t_INT: return mulsi(s,x);
2853 70335314 : case t_REAL: return s? mulsr(s,x): gen_0; /* gmul semantic */
2854 174365 : case t_INTMOD: { GEN p = gel(x,1);
2855 174365 : z = cgetg(3,t_INTMOD);
2856 174366 : gel(z,2) = gc_INT((pari_sp)z, modii(mulsi(s,gel(x,2)), p));
2857 174365 : gel(z,1) = icopy(p); return z;
2858 : }
2859 548098 : case t_FFELT: return FF_Z_mul(x,stoi(s));
2860 6217613 : case t_FRAC:
2861 6217613 : if (!s) return gen_0;
2862 6202241 : z = cgetg(3,t_FRAC);
2863 6202241 : i = labs(s); i = ugcdiu(gel(x,2), i);
2864 6202241 : if (i == 1)
2865 : {
2866 2270597 : gel(z,2) = icopy(gel(x,2));
2867 2270597 : gel(z,1) = mulis(gel(x,1), s);
2868 : }
2869 : else
2870 : {
2871 3931644 : gel(z,2) = diviuexact(gel(x,2), (ulong)i);
2872 3931644 : gel(z,1) = mulis(gel(x,1), s/i);
2873 3931644 : fix_frac_if_int(z);
2874 : }
2875 6202241 : return z;
2876 :
2877 14522485 : case t_COMPLEX:
2878 14522485 : if (!s) return gen_0;
2879 13396681 : z = cgetg(3, t_COMPLEX);
2880 13396674 : gel(z,1) = gmulsg(s,gel(x,1));
2881 13394979 : gel(z,2) = gmulsg(s,gel(x,2)); return z;
2882 :
2883 1449 : case t_PADIC:
2884 1449 : if (!s) return gen_0;
2885 1449 : av = avma; return gc_upto(av, mulpp(cvtop2(stoi(s),x), x));
2886 :
2887 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2888 7 : gel(z,1) = ZX_copy(gel(x,1));
2889 7 : gel(z,2) = gmulsg(s,gel(x,2));
2890 7 : gel(z,3) = gmulsg(s,gel(x,3)); return z;
2891 :
2892 208516 : case t_POLMOD:
2893 208516 : retmkpolmod(gmulsg(s,gel(x,2)), RgX_copy(gel(x,1)));
2894 :
2895 821368 : case t_POL:
2896 821368 : if (!signe(x)) return RgX_copy(x);
2897 800991 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2898 3327613 : pari_APPLY_pol(gmulsg(s,gel(x,i)));
2899 :
2900 182 : case t_SER:
2901 182 : if (ser_isexactzero(x)) return gcopy(x);
2902 182 : if (!s) return Rg_get_0(x);
2903 3864 : pari_APPLY_ser(gmulsg(s,gel(x,i)));
2904 :
2905 0 : case t_RFRAC:
2906 0 : if (!s) return zeropol(varn(gel(x,2)));
2907 0 : if (s == 1) return gcopy(x);
2908 0 : if (s == -1) return gneg(x);
2909 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), stoi(s));
2910 :
2911 1304584 : case t_VEC: case t_COL: case t_MAT:
2912 4109862 : pari_APPLY_same(gmulsg(s,gel(x,i)));
2913 : }
2914 0 : pari_err_TYPE("gmulsg",x);
2915 : return NULL; /* LCOV_EXCL_LINE */
2916 : }
2917 :
2918 : GEN
2919 124385044 : gmulug(ulong s, GEN x)
2920 : {
2921 : pari_sp av;
2922 : long i;
2923 : GEN z;
2924 :
2925 124385044 : switch(typ(x))
2926 : {
2927 122251738 : case t_INT: return mului(s,x);
2928 1064579 : case t_REAL: return s? mulur(s,x): gen_0; /* gmul semantic */
2929 364 : case t_INTMOD: { GEN p = gel(x,1);
2930 364 : z = cgetg(3,t_INTMOD);
2931 364 : gel(z,2) = gc_INT((pari_sp)z, modii(mului(s,gel(x,2)), p));
2932 364 : gel(z,1) = icopy(p); return z;
2933 : }
2934 413 : case t_FFELT: return FF_Z_mul(x,utoi(s));
2935 981585 : case t_FRAC:
2936 981585 : if (!s) return gen_0;
2937 981571 : z = cgetg(3,t_FRAC);
2938 981571 : i = ugcdiu(gel(x,2), s);
2939 981571 : if (i == 1)
2940 : {
2941 822559 : gel(z,2) = icopy(gel(x,2));
2942 822559 : gel(z,1) = muliu(gel(x,1), s);
2943 : }
2944 : else
2945 : {
2946 159012 : gel(z,2) = diviuexact(gel(x,2), i);
2947 159012 : gel(z,1) = muliu(gel(x,1), s/i);
2948 159012 : fix_frac_if_int(z);
2949 : }
2950 981571 : return z;
2951 :
2952 35945 : case t_COMPLEX:
2953 35945 : if (!s) return gen_0;
2954 35945 : z = cgetg(3, t_COMPLEX);
2955 35945 : gel(z,1) = gmulug(s,gel(x,1));
2956 35945 : gel(z,2) = gmulug(s,gel(x,2)); return z;
2957 :
2958 7342 : case t_PADIC:
2959 7342 : if (!s) return gen_0;
2960 7342 : av = avma; return gc_upto(av, mulpp(cvtop2(utoi(s),x), x));
2961 :
2962 0 : case t_QUAD: z = cgetg(4, t_QUAD);
2963 0 : gel(z,1) = ZX_copy(gel(x,1));
2964 0 : gel(z,2) = gmulug(s,gel(x,2));
2965 0 : gel(z,3) = gmulug(s,gel(x,3)); return z;
2966 :
2967 6363 : case t_POLMOD:
2968 6363 : retmkpolmod(gmulug(s,gel(x,2)), RgX_copy(gel(x,1)));
2969 :
2970 19096 : case t_POL:
2971 19096 : if (!signe(x)) return RgX_copy(x);
2972 18319 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2973 61873 : pari_APPLY_pol(gmulug(s,gel(x,i)));
2974 :
2975 0 : case t_SER:
2976 0 : if (ser_isexactzero(x)) return gcopy(x);
2977 0 : if (!s) return Rg_get_0(x);
2978 0 : pari_APPLY_ser(gmulug(s,gel(x,i)));
2979 :
2980 0 : case t_RFRAC:
2981 0 : if (!s) return zeropol(varn(gel(x,2)));
2982 0 : if (s == 1) return gcopy(x);
2983 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), utoi(s));
2984 :
2985 17619 : case t_VEC: case t_COL: case t_MAT:
2986 96495 : pari_APPLY_same(gmulug(s,gel(x,i)));
2987 : }
2988 0 : pari_err_TYPE("gmulsg",x);
2989 : return NULL; /* LCOV_EXCL_LINE */
2990 : }
2991 :
2992 : /********************************************************************/
2993 : /** **/
2994 : /** SIMPLE DIVISION **/
2995 : /** **/
2996 : /********************************************************************/
2997 :
2998 : GEN
2999 13123128 : gdivgs(GEN x, long s)
3000 : {
3001 13123128 : long tx = typ(x), i;
3002 : pari_sp av;
3003 : GEN z;
3004 :
3005 13123128 : if (!s)
3006 : {
3007 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3008 0 : pari_err_INV("gdivgs",gen_0);
3009 : }
3010 13123130 : switch(tx)
3011 : {
3012 1591257 : case t_INT: return Qdivis(x, s);
3013 8649989 : case t_REAL: return divrs(x,s);
3014 :
3015 357 : case t_INTMOD:
3016 357 : z = cgetg(3, t_INTMOD);
3017 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3018 :
3019 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3020 :
3021 556780 : case t_FRAC: z = cgetg(3, t_FRAC);
3022 556780 : i = ugcdiu(gel(x,1), labs(s));
3023 556780 : if (i == 1)
3024 : {
3025 413362 : gel(z,2) = mulsi(s, gel(x,2));
3026 413362 : gel(z,1) = icopy(gel(x,1));
3027 : }
3028 : else
3029 : {
3030 143418 : gel(z,2) = mulsi(s/i, gel(x,2));
3031 143418 : gel(z,1) = divis(gel(x,1), i);
3032 : }
3033 556780 : normalize_frac(z);
3034 556780 : fix_frac_if_int(z); return z;
3035 :
3036 1841806 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3037 1841806 : gel(z,1) = gdivgs(gel(x,1),s);
3038 1841805 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3039 :
3040 133 : case t_PADIC: /* divpT */
3041 : {
3042 133 : GEN p = padic_p(x);
3043 133 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3044 133 : av = avma;
3045 133 : return gc_upto(av, divpp(x, cvtop2(stoi(s),x)));
3046 : }
3047 :
3048 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3049 28 : gel(z,1) = ZX_copy(gel(x,1));
3050 28 : gel(z,2) = gdivgs(gel(x,2),s);
3051 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3052 :
3053 37282 : case t_POLMOD:
3054 37282 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3055 :
3056 91 : case t_RFRAC:
3057 91 : if (s == 1) return gcopy(x);
3058 84 : else if (s == -1) return gneg(x);
3059 84 : return div_rfrac_scal(x, stoi(s));
3060 :
3061 156466 : case t_POL: pari_APPLY_pol_normalized(gdivgs(gel(x,i),s));
3062 0 : case t_SER: pari_APPLY_ser_normalized(gdivgs(gel(x,i),s));
3063 407593 : case t_VEC:
3064 : case t_COL:
3065 894977 : case t_MAT: pari_APPLY_same(gdivgs(gel(x,i),s));
3066 : }
3067 0 : pari_err_TYPE2("/",x, stoi(s));
3068 : return NULL; /* LCOV_EXCL_LINE */
3069 : }
3070 :
3071 : GEN
3072 54785016 : gdivgu(GEN x, ulong s)
3073 : {
3074 54785016 : long tx = typ(x), i;
3075 : pari_sp av;
3076 : GEN z;
3077 :
3078 54785016 : if (!s)
3079 : {
3080 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3081 0 : pari_err_INV("gdivgu",gen_0);
3082 : }
3083 54784913 : switch(tx)
3084 : {
3085 18077323 : case t_INT: return Qdiviu(x, s);
3086 12895672 : case t_REAL: return divru(x,s);
3087 :
3088 210315 : case t_INTMOD:
3089 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3090 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3091 :
3092 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3093 :
3094 643003 : case t_FRAC: z = cgetg(3, t_FRAC);
3095 643003 : i = ugcdiu(gel(x,1), s);
3096 643003 : if (i == 1)
3097 : {
3098 463131 : gel(z,2) = mului(s, gel(x,2));
3099 463131 : gel(z,1) = icopy(gel(x,1));
3100 : }
3101 : else
3102 : {
3103 179872 : gel(z,2) = mului(s/i, gel(x,2));
3104 179872 : gel(z,1) = divis(gel(x,1), i);
3105 : }
3106 643003 : normalize_frac(z);
3107 643003 : fix_frac_if_int(z); return z;
3108 :
3109 11347245 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3110 11347245 : gel(z,1) = gdivgu(gel(x,1),s);
3111 11347246 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3112 :
3113 11551708 : case t_PADIC: /* divpT */
3114 : {
3115 11551708 : GEN p = padic_p(x);
3116 11551708 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3117 11416068 : av = avma;
3118 11416068 : return gc_upto(av, divpp(x, cvtop2(utoi(s),x)));
3119 : }
3120 :
3121 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3122 0 : gel(z,1) = ZX_copy(gel(x,1));
3123 0 : gel(z,2) = gdivgu(gel(x,2),s);
3124 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3125 :
3126 1456 : case t_POLMOD:
3127 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3128 :
3129 56 : case t_RFRAC:
3130 56 : if (s == 1) return gcopy(x);
3131 56 : return div_rfrac_scal(x, utoi(s));
3132 :
3133 329287 : case t_POL: pari_APPLY_pol_normalized(gdivgu(gel(x,i),s));
3134 16786 : case t_SER: pari_APPLY_ser_normalized(gdivgu(gel(x,i),s));
3135 336 : case t_VEC:
3136 : case t_COL:
3137 1162 : case t_MAT: pari_APPLY_same(gdivgu(gel(x,i),s));
3138 : }
3139 0 : pari_err_TYPE2("/",x, utoi(s));
3140 : return NULL; /* LCOV_EXCL_LINE */
3141 : }
3142 :
3143 : /* x / (i*(i+1)) */
3144 : GEN
3145 224483219 : divrunextu(GEN x, ulong i)
3146 : {
3147 224483219 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3148 0 : return divri(x, muluu(i , i+1));
3149 : else
3150 224483219 : return divru(x, i*(i+1));
3151 : }
3152 : /* x / (i*(i+1)) */
3153 : GEN
3154 1583873 : gdivgunextu(GEN x, ulong i)
3155 : {
3156 1583873 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3157 0 : return gdivgu(x, i*(i+1));
3158 : else
3159 1583873 : return gdiv(x, muluu(i, i+1));
3160 : }
3161 :
3162 : /* True shift (exact multiplication by 2^n) */
3163 : GEN
3164 138467601 : gmul2n(GEN x, long n)
3165 : {
3166 : GEN z, a, b;
3167 : long k, l;
3168 :
3169 138467601 : switch(typ(x))
3170 : {
3171 38937884 : case t_INT:
3172 38937884 : if (n>=0) return shifti(x,n);
3173 7451361 : if (!signe(x)) return gen_0;
3174 4887256 : l = vali(x); n = -n;
3175 4887296 : if (n<=l) return shifti(x,-n);
3176 412731 : z = cgetg(3,t_FRAC);
3177 412731 : gel(z,1) = shifti(x,-l);
3178 412731 : gel(z,2) = int2n(n-l); return z;
3179 :
3180 63442428 : case t_REAL:
3181 63442428 : return shiftr(x,n);
3182 :
3183 180965 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3184 180965 : z = cgetg(3,t_INTMOD);
3185 180965 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3186 82804 : gel(z,2) = gc_INT((pari_sp)z, modii(shifti(a,n), b));
3187 82802 : gel(z,1) = icopy(b); return z;
3188 :
3189 217424 : case t_FFELT: return FF_mul2n(x,n);
3190 :
3191 1309851 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3192 1309851 : l = vali(a);
3193 1309851 : k = vali(b);
3194 1309851 : if (n+l >= k)
3195 : {
3196 425494 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3197 331998 : l = n-k; k = -k;
3198 : }
3199 : else
3200 : {
3201 884357 : k = -(l+n); l = -l;
3202 : }
3203 1216355 : z = cgetg(3,t_FRAC);
3204 1216355 : gel(z,1) = shifti(a,l);
3205 1216355 : gel(z,2) = shifti(b,k); return z;
3206 :
3207 10765897 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3208 10765893 : gel(z,1) = gmul2n(gel(x,1),n);
3209 10765898 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3210 :
3211 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3212 105 : gel(z,1) = ZX_copy(gel(x,1));
3213 105 : gel(z,2) = gmul2n(gel(x,2),n);
3214 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3215 :
3216 177009 : case t_POLMOD:
3217 177009 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3218 :
3219 6750788 : case t_POL:
3220 24029040 : pari_APPLY_pol(gmul2n(gel(x,i),n));
3221 102875 : case t_SER:
3222 102875 : if (ser_isexactzero(x)) return gcopy(x);
3223 634312 : pari_APPLY_ser(gmul2n(gel(x,i),n));
3224 16576653 : case t_VEC: case t_COL: case t_MAT:
3225 60986455 : pari_APPLY_same(gmul2n(gel(x,i),n));
3226 :
3227 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3228 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3229 :
3230 7644 : case t_PADIC: /* int2n wrong if n < 0 */
3231 7644 : return gmul(gmul2n(gen_1,n),x);
3232 : }
3233 0 : pari_err_TYPE("gmul2n",x);
3234 : return NULL; /* LCOV_EXCL_LINE */
3235 : }
3236 :
3237 : /*******************************************************************/
3238 : /* */
3239 : /* INVERSE */
3240 : /* */
3241 : /*******************************************************************/
3242 : static GEN
3243 189625 : inv_polmod(GEN T, GEN x)
3244 : {
3245 189625 : GEN z = cgetg(3,t_POLMOD), a;
3246 189624 : gel(z,1) = RgX_copy(T);
3247 189622 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3248 64912 : a = ginv(x);
3249 : else
3250 : {
3251 124710 : if (lg(T) == 5) /* quadratic fields */
3252 13125 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3253 : else
3254 111585 : a = RgXQ_inv(x, T);
3255 : }
3256 189622 : gel(z,2) = a; return z;
3257 : }
3258 :
3259 : GEN
3260 36577061 : ginv(GEN x)
3261 : {
3262 : long s;
3263 : pari_sp av;
3264 : GEN z, y;
3265 :
3266 36577061 : switch(typ(x))
3267 : {
3268 9508099 : case t_INT:
3269 9508099 : if (is_pm1(x)) return icopy(x);
3270 7958581 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3271 7958567 : z = cgetg(3,t_FRAC);
3272 7964847 : gel(z,1) = s<0? gen_m1: gen_1;
3273 7964847 : gel(z,2) = absi(x); return z;
3274 :
3275 10862182 : case t_REAL: return invr(x);
3276 :
3277 5733 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3278 5733 : gel(z,1) = icopy(gel(x,1));
3279 5733 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3280 :
3281 532326 : case t_FRAC: {
3282 532326 : GEN a = gel(x,1), b = gel(x,2);
3283 532326 : s = signe(a);
3284 532326 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3285 247497 : z = cgetg(3,t_FRAC);
3286 247497 : gel(z,1) = icopy(b);
3287 247497 : gel(z,2) = icopy(a);
3288 247497 : normalize_frac(z); return z;
3289 : }
3290 9618536 : case t_COMPLEX:
3291 9618536 : av = avma;
3292 9618536 : return gc_upto(av, divcR(conj_i(x), cxnorm(x)));
3293 :
3294 273 : case t_QUAD:
3295 273 : av = avma;
3296 273 : return gc_upto(av, gdiv(conj_i(x), quadnorm(x)));
3297 :
3298 358358 : case t_PADIC:
3299 : {
3300 358358 : GEN p = padic_p(x), pd = padic_pd(x), u = padic_u(x);
3301 358358 : long d = precp(x);
3302 358358 : if (!signe(u)) pari_err_INV("ginv",x);
3303 358351 : retmkpadic(Zp_inv(u, p, d), icopy(p), icopy(pd), -valp(x), d);
3304 : }
3305 :
3306 189625 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3307 14060 : case t_FFELT: return FF_inv(x);
3308 5322037 : case t_POL: return gred_rfrac_simple(gen_1,x);
3309 26369 : case t_SER: return ser_inv(x);
3310 2933 : case t_RFRAC:
3311 : {
3312 2933 : GEN n = gel(x,1), d = gel(x,2);
3313 2933 : pari_sp av = avma, ltop;
3314 2933 : if (gequal0(n)) pari_err_INV("ginv",x);
3315 :
3316 2933 : n = simplify_shallow(n);
3317 2933 : if (typ(n) != t_POL || varn(n) != varn(d))
3318 : {
3319 2933 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3320 679 : ltop = avma;
3321 679 : z = RgX_Rg_div(d,n);
3322 : } else {
3323 0 : ltop = avma;
3324 0 : z = cgetg(3,t_RFRAC);
3325 0 : gel(z,1) = RgX_copy(d);
3326 0 : gel(z,2) = RgX_copy(n);
3327 : }
3328 679 : stackdummy(av, ltop);
3329 679 : return z;
3330 : }
3331 :
3332 21 : case t_VEC: if (!is_ext_qfr(x)) break;
3333 : case t_QFB:
3334 99002 : return qfbpow(x, gen_m1);
3335 38939 : case t_MAT:
3336 38939 : y = RgM_inv(x);
3337 38932 : if (!y) pari_err_INV("ginv",x);
3338 38862 : return y;
3339 28 : case t_VECSMALL:
3340 : {
3341 28 : long i, lx = lg(x)-1;
3342 28 : y = zero_zv(lx);
3343 112 : for (i=1; i<=lx; i++)
3344 : {
3345 84 : long xi = x[i];
3346 84 : if (xi<1 || xi>lx || y[xi])
3347 0 : pari_err_TYPE("ginv [not a permutation]", x);
3348 84 : y[xi] = i;
3349 : }
3350 28 : return y;
3351 : }
3352 : }
3353 6 : pari_err_TYPE("inverse",x);
3354 : return NULL; /* LCOV_EXCL_LINE */
3355 : }
|
