Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** GENERIC OPERATIONS **/
18 : /** (first part) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_mod
25 :
26 : /* assume z[1] was created last */
27 : #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
28 : z = gc_upto((pari_sp)(z+3), gel(z,1));
29 :
30 : /* assume z[1] was created last */
31 : #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
32 : z = gc_upto((pari_sp)(z+3), gel(z,1));\
33 : else\
34 : gc_slice_unsafe((pari_sp)z, tetpil, z+1, 2); }
35 :
36 : static void
37 112 : warn_coercion(GEN x, GEN y, GEN z)
38 : {
39 112 : if (DEBUGLEVEL)
40 56 : pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
41 112 : }
42 :
43 : static long
44 35 : kro_quad(GEN x, GEN y)
45 35 : { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
46 :
47 : /* is -1 not a square in Zp, assume p prime */
48 : INLINE int
49 42 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
50 :
51 : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
52 : static GEN mulpp(GEN x, GEN y);
53 : static GEN divpp(GEN x, GEN y);
54 : /* Argument codes for inline routines
55 : * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
56 : * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
57 : * T: some type (to be converted to PADIC)
58 : */
59 : static GEN
60 313029612 : addRc(GEN x, GEN y) {
61 313029612 : GEN z = cgetg(3,t_COMPLEX);
62 313013960 : gel(z,1) = gadd(x,gel(y,1));
63 313015470 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 358400079 : mulRc(GEN x, GEN y)
67 : {
68 358400079 : GEN a = gel(y,1), b = gel(y,2), z = cgetg(3,t_COMPLEX);
69 358438830 : gel(z,1) = (typ(x) != t_INTMOD && isintzero(a))? gen_0: gmul(x, a);
70 358671034 : gel(z,2) = gmul(x, b); return z;
71 : }
72 : static GEN
73 3408011 : divRc(GEN x, GEN y)
74 : {
75 3408011 : GEN a = gel(y,1), b = gel(y,2), z = cgetg(3,t_COMPLEX);
76 3408009 : pari_sp av = avma, av2;
77 3408009 : if (isintzero(a) && typ(x) != t_INTMOD)
78 : {
79 526099 : gel(z,1) = gen_0;
80 526099 : gel(z,2) = gc_upto(av, gdiv(x, gneg(b)));
81 : }
82 : else
83 : {
84 2881910 : GEN t = gdiv(x, cxnorm(y)), mb = gneg(b);
85 2881908 : av2 = avma;
86 2881908 : gel(z,1) = gmul(t, a);
87 2881899 : gel(z,2) = gmul(t, mb);
88 2881894 : gc_slice_unsafe(av, av2, z+1, 2);
89 : }
90 3408004 : 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 24801170 : divcR(GEN x, GEN y) {
100 24801170 : GEN z = cgetg(3,t_COMPLEX);
101 24801044 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
102 24800842 : gel(z,2) = gdiv(gel(x,2), y); return z;
103 : }
104 : static GEN
105 1596 : addRq(GEN x, GEN y) {
106 1596 : GEN z = cgetg(4,t_QUAD);
107 1596 : gel(z,1) = ZX_copy(gel(y,1));
108 1596 : gel(z,2) = gadd(x, gel(y,2));
109 1596 : 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 38167471 : mulrfrac(GEN x, GEN y)
126 : {
127 : pari_sp av;
128 38167471 : GEN z, a = gel(y,1), b = gel(y,2);
129 38167471 : if (is_pm1(a)) /* frequent special case */
130 : {
131 12826446 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
132 12823644 : return z;
133 : }
134 25340753 : av = avma;
135 25340753 : 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 4985910 : mulTp(GEN x, GEN y) { pari_sp av = avma;
163 4985910 : return gc_upto(av, mulpp(cvtop2(x,y), y));
164 : }
165 : /* y PADIC, non zero x / y by converting x to padic */
166 : static GEN
167 3674 : divTp(GEN x, GEN y) { pari_sp av = avma;
168 3674 : return gc_upto(av, divpp(cvtop2(x,y), y));
169 : }
170 : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
171 : * converted to O(p^e) and division by 0 */
172 : static GEN
173 1204910 : divpT(GEN x, GEN y) { pari_sp av = avma;
174 1204910 : return gc_upto(av, divpp(x, cvtop2(y,x)));
175 : }
176 :
177 : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
178 : * clean memory from z on */
179 : static GEN
180 3416149 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 3416149 : if (lgefint(X) == 3) {
182 3091895 : ulong u = Fl_add(itou(x),itou(y), X[2]);
183 3091895 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
184 : }
185 : else {
186 324254 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
187 324251 : gel(z,2) = gc_INT((pari_sp)z, u);
188 : }
189 3416148 : gel(z,1) = icopy(X); return z;
190 : }
191 : static GEN
192 1031401 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
193 1031401 : if (lgefint(X) == 3) {
194 658949 : ulong u = Fl_sub(itou(x),itou(y), X[2]);
195 658949 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
196 : }
197 : else {
198 372452 : GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
199 372452 : gel(z,2) = gc_INT((pari_sp)z, u);
200 : }
201 1031401 : gel(z,1) = icopy(X); return z;
202 : }
203 : /* cf add_intmod_same */
204 : static GEN
205 3138696 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
206 3138696 : if (lgefint(X) == 3) {
207 2996214 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
208 2996214 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 142482 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x,y), X) );
212 3138700 : gel(z,1) = icopy(X); return z;
213 : }
214 : /* cf add_intmod_same */
215 : static GEN
216 341500 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
217 : {
218 341500 : if (lgefint(X) == 3) {
219 318971 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
220 318957 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
221 : }
222 : else
223 22529 : gel(z,2) = gc_INT((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
224 341486 : 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 9959295 : rfrac_denom_mul_scal(GEN d, GEN y)
237 : {
238 9959295 : GEN D = RgX_Rg_mul(d, y);
239 9959295 : 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 9959295 : return D;
245 : }
246 :
247 : static int
248 57792183 : Leading_is_neg(GEN x)
249 : {
250 122218150 : while (typ(x) == t_POL) x = leading_coeff(x);
251 57792183 : return is_real_t(typ(x))? gsigne(x) < 0: 0;
252 : }
253 :
254 : /* x = Rg_get_1(something); shall we use it for base change ?
255 : * Not if gen_1 (over Z) or t_PADIC (over Zp, padic prec is independent in
256 : * each coefficient) */
257 : static int
258 154407967 : base_change(GEN x) { return x != gen_1 && typ(x) != t_PADIC; }
259 :
260 : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
261 : GEN
262 57797147 : gred_rfrac_simple(GEN n, GEN d)
263 : {
264 : GEN _1n, _1d, c, cn, cd, z;
265 57797147 : long dd = degpol(d);
266 :
267 57797147 : if (dd <= 0)
268 : {
269 4957 : if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
270 4957 : n = gdiv(n, gel(d,2));
271 4957 : if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
272 4957 : return n;
273 : }
274 57792190 : _1n = Rg_get_1(n);
275 57792190 : _1d = Rg_get_1(d);
276 57792190 : if (base_change(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
277 57792190 : if (base_change(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
278 57792183 : if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
279 57792183 : cd = content(d);
280 59654048 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
281 57792183 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
282 57792183 : if (!gequal1(cd)) {
283 6586902 : d = RgX_Rg_div(d,cd);
284 6586902 : if (!gequal1(cn))
285 : {
286 1313733 : if (gequal0(cn)) {
287 91 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
288 0 : n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
289 0 : c = gen_1;
290 : } else {
291 1313642 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
292 1313642 : c = gdiv(cn,cd);
293 : }
294 : }
295 : else
296 5273169 : c = ginv(cd);
297 : } else {
298 51205281 : if (!gequal1(cn))
299 : {
300 3341290 : if (gequal0(cn)) {
301 966 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
302 21 : c = gen_1;
303 : } else {
304 3340324 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
305 3340324 : c = cn;
306 : }
307 : } else {
308 47863991 : GEN y = cgetg(3,t_RFRAC);
309 47863991 : gel(y,1) = gcopy(n);
310 47863991 : gel(y,2) = RgX_copy(d); return y;
311 : }
312 : }
313 :
314 9927156 : if (typ(c) == t_POL)
315 : {
316 912651 : z = c;
317 952698 : do { z = content(z); } while (typ(z) == t_POL);
318 912651 : cd = denom_i(z);
319 912651 : cn = gmul(c, cd);
320 : }
321 : else
322 : {
323 9014505 : cn = numer_i(c);
324 9014505 : cd = denom_i(c);
325 : }
326 9927156 : z = cgetg(3,t_RFRAC);
327 9927156 : gel(z,1) = gmul(n, cn);
328 9927156 : gel(z,2) = d = rfrac_denom_mul_scal(d, cd);
329 : /* can happen: Pol(O(17^-1)) / Pol([Mod(9,23), O(23^-3)]) */
330 9927156 : if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
331 9927156 : return z;
332 : }
333 :
334 : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
335 : static GEN
336 207018 : fix_rfrac(GEN x, long d)
337 : {
338 : GEN z, N, D;
339 207018 : if (!d || typ(x) == t_POL) return x;
340 165697 : z = cgetg(3, t_RFRAC);
341 165697 : N = gel(x,1);
342 165697 : D = gel(x,2);
343 165697 : if (d > 0) {
344 6650 : gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
345 160944 : : monomialcopy(N,d,varn(D));
346 160881 : gel(z, 2) = RgX_copy(D);
347 : } else {
348 4816 : gel(z, 1) = gcopy(N);
349 4816 : gel(z, 2) = RgX_shift(D, -d);
350 : }
351 165697 : return z;
352 : }
353 :
354 : /* d a t_POL; simplify the rational function n / d */
355 : static GEN
356 44825263 : gred_rfrac(GEN n, GEN d)
357 : {
358 : GEN y, z, _1n, _1d;
359 : long v, vd, vn;
360 :
361 44825263 : n = simplify_shallow(n);
362 44825263 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
363 37968185 : d = simplify_shallow(d);
364 37968185 : if (typ(d) != t_POL) return gdiv(n,d);
365 36791289 : vd = varn(d);
366 36791289 : if (typ(n) != t_POL)
367 : {
368 20654411 : if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
369 20653004 : if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
370 0 : pari_err_BUG("gred_rfrac [incompatible variables]");
371 : }
372 16136878 : vn = varn(n);
373 16136878 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
374 15995081 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
375 15835308 : _1n = Rg_get_1(n);
376 15835308 : _1d = Rg_get_1(d);
377 15835308 : if (base_change(_1n) && !gidentical(_1n, _1d))
378 : {
379 168 : d = gmul(d, _1n);
380 168 : if (!signe(d)) pari_err_INV("gdiv",d);
381 : }
382 15835308 : if (base_change(_1d) && !gidentical(_1n, _1d))
383 : {
384 490 : n = gmul(n, _1d);
385 490 : if (!signe(n)) return gdiv(simplify_shallow(n),d);
386 : }
387 : /* now n and d are t_POLs in the same variable */
388 15835294 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
389 15835294 : if (!degpol(d))
390 : {
391 12798432 : n = RgX_Rg_div(n,gel(d,2));
392 12798432 : return v? RgX_mulXn(n,v): n;
393 : }
394 :
395 : /* X does not divide gcd(n,d), deg(d) > 0 */
396 3036862 : if (!isinexact(n) && !isinexact(d))
397 : {
398 3036617 : y = RgX_divrem(n, d, &z);
399 3036617 : if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
400 206773 : z = RgX_gcd(d, z);
401 206773 : if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
402 : }
403 207018 : return fix_rfrac(gred_rfrac_simple(n,d), v);
404 : }
405 :
406 : /* x,y t_INT, return x/y in reduced form */
407 : GEN
408 144868112 : Qdivii(GEN x, GEN y)
409 : {
410 144868112 : pari_sp av = avma;
411 : GEN r, q;
412 :
413 144868112 : if (lgefint(y) == 3)
414 : {
415 124813987 : q = Qdiviu(x, y[2]);
416 124813594 : if (signe(y) > 0) return q;
417 11101416 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
418 11101431 : return q;
419 : }
420 20054125 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
421 20054878 : if (equali1(x))
422 : {
423 5186801 : if (!signe(y)) pari_err_INV("gdiv",y);
424 5186689 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
425 : }
426 14867889 : q = dvmdii(x,y,&r);
427 14868261 : if (r == gen_0) return q; /* gen_0 intended */
428 8930090 : r = gcdii(y, r);
429 8930000 : if (lgefint(r) == 3)
430 : {
431 8015685 : ulong t = r[2];
432 8015685 : set_avma(av);
433 8015611 : if (t == 1) q = mkfraccopy(x,y);
434 : else
435 : {
436 2886265 : q = cgetg(3,t_FRAC);
437 2886619 : gel(q,1) = diviuexact(x,t);
438 2886537 : gel(q,2) = diviuexact(y,t);
439 : }
440 : }
441 : else
442 : { /* rare: r and q left on stack for efficiency */
443 914315 : q = cgetg(3,t_FRAC);
444 914323 : gel(q,1) = diviiexact(x,r);
445 914323 : gel(q,2) = diviiexact(y,r);
446 : }
447 8930198 : normalize_frac(q); return q;
448 : }
449 :
450 : static GEN
451 65599229 : utoi_sign(ulong x, long s) { return s > 0? utoipos(x): utoineg(x); }
452 :
453 : /* x t_INT, return x/y in reduced form */
454 : GEN
455 152508083 : Qdiviu(GEN x, ulong y)
456 : {
457 : pari_sp av;
458 : ulong r, t;
459 : long s;
460 : GEN q;
461 :
462 152508083 : if (y == 1) return icopy(x);
463 130812754 : if (!y) pari_err_INV("gdiv",gen_0);
464 130815873 : s = signe(x); if (!s) return gen_0;
465 83537070 : if (lgefint(x) == 3)
466 : {
467 77441794 : ulong xx = uel(x,2), g;
468 77441794 : if (xx == 1) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
469 65599376 : g = ugcd(xx, y);
470 65598846 : if (g != 1)
471 : {
472 43698971 : xx /= g; if (y == g) return utoi_sign(xx, s);
473 17231401 : y /= g;
474 : }
475 39131276 : retmkfrac(utoi_sign(xx, s), utoipos(y));
476 : }
477 6095276 : av = avma; q = absdiviu_rem(x,y,&r);
478 6095321 : if (!r)
479 : {
480 2044310 : if (s < 0) togglesign(q);
481 2044295 : return q;
482 : }
483 4051011 : t = ugcd(y, r); set_avma(av);
484 4051022 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
485 1486051 : retmkfrac(diviuexact(x,t), utoipos(y / t));
486 : }
487 :
488 : /* x t_INT, return x/y in reduced form */
489 : GEN
490 2582080 : Qdivis(GEN x, long y)
491 : {
492 : GEN q;
493 2582080 : if (y >= 0) return Qdiviu(x, y);
494 181428 : q = Qdiviu(x, -y);
495 181432 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
496 181432 : return q;
497 : }
498 :
499 : /*******************************************************************/
500 : /* */
501 : /* CONJUGATION */
502 : /* */
503 : /*******************************************************************/
504 : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
505 : static GEN
506 18725 : quad_polmod_conj(GEN x, GEN y)
507 : {
508 : GEN z, u, v, a, b;
509 18725 : if (typ(x) != t_POL) return gcopy(x);
510 18725 : if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
511 18725 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
512 18725 : b = gel(y,3); v = gel(x,2);
513 18725 : z = cgetg(4, t_POL); z[1] = x[1];
514 18725 : gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
515 18725 : gel(z,3) = gneg(u); return z;
516 : }
517 : static GEN
518 18725 : quad_polmod_norm(GEN x, GEN y)
519 : {
520 : GEN z, u, v, a, b, c;
521 18725 : if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
522 18725 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
523 18725 : b = gel(y,3); v = gel(x,2);
524 18725 : c = gel(y,2);
525 18725 : z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
526 18725 : if (!gequal1(a)) z = gdiv(z, a);
527 18725 : return gadd(z, gsqr(v));
528 : }
529 :
530 : GEN
531 32592171 : conj_i(GEN x)
532 : {
533 32592171 : switch(typ(x))
534 : {
535 7465710 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
536 7465710 : return x;
537 :
538 24961218 : case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
539 :
540 1190 : case t_QUAD:
541 : {
542 1190 : GEN y = cgetg(4,t_QUAD);
543 1190 : gel(y,1) = gel(x,1);
544 1190 : gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
545 1190 : : gadd(gel(x,2), gel(x,3));
546 1190 : gel(y,3) = gneg(gel(x,3)); return y;
547 : }
548 350 : case t_POL: pari_APPLY_pol_normalized(conj_i(gel(x,i)));
549 33047 : case t_SER: pari_APPLY_ser_normalized(conj_i(gel(x,i)));
550 :
551 153601 : case t_RFRAC:
552 : case t_VEC:
553 : case t_COL:
554 555555 : case t_MAT: pari_APPLY_same(conj_i(gel(x,i)));
555 :
556 0 : case t_POLMOD:
557 : {
558 0 : GEN X = gel(x,1);
559 0 : long d = degpol(X);
560 0 : if (d < 2) return x;
561 0 : if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
562 : }
563 : }
564 0 : pari_err_TYPE("gconj",x);
565 : return NULL; /* LCOV_EXCL_LINE */
566 : }
567 : GEN
568 390611 : gconj(GEN x)
569 390611 : { pari_sp av = avma; return gc_GEN(av, conj_i(x)); }
570 :
571 : GEN
572 84 : conjvec(GEN x,long prec)
573 : {
574 : long lx, s, i;
575 : GEN z;
576 :
577 84 : switch(typ(x))
578 : {
579 0 : case t_INT: case t_INTMOD: case t_FRAC:
580 0 : return mkcolcopy(x);
581 :
582 0 : case t_COMPLEX: case t_QUAD:
583 0 : z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
584 :
585 28 : case t_FFELT:
586 28 : return FF_conjvec(x);
587 :
588 0 : case t_VEC: case t_COL:
589 0 : lx = lg(x); z = cgetg(lx,t_MAT);
590 0 : if (lx == 1) return z;
591 0 : gel(z,1) = conjvec(gel(x,1),prec);
592 0 : s = lgcols(z);
593 0 : for (i=2; i<lx; i++)
594 : {
595 0 : gel(z,i) = conjvec(gel(x,i),prec);
596 0 : if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
597 : }
598 0 : break;
599 :
600 56 : case t_POLMOD: {
601 56 : GEN T = gel(x,1), r;
602 : pari_sp av;
603 :
604 56 : lx = lg(T);
605 56 : if (lx <= 3) return cgetg(1,t_COL);
606 56 : x = gel(x,2);
607 238 : for (i=2; i<lx; i++)
608 : {
609 189 : GEN c = gel(T,i);
610 189 : switch(typ(c)) {
611 7 : case t_INTMOD: {
612 7 : GEN p = gel(c,1);
613 : pari_sp av;
614 7 : if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
615 7 : av = avma;
616 7 : T = RgX_to_FpX(T,p);
617 7 : x = RgX_to_FpX(x, p);
618 7 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
619 7 : z = FpXQC_to_mod(FpXQ_conjvec(x, T, p), T, p);
620 7 : return gc_upto(av, z);
621 : }
622 182 : case t_INT:
623 182 : case t_FRAC: break;
624 0 : default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
625 : }
626 : }
627 49 : if (typ(x) != t_POL)
628 : {
629 0 : if (!is_rational_t(typ(x)))
630 0 : pari_err_TYPE("conjvec [not a rational t_POL]",x);
631 0 : retconst_col(lx-3, gcopy(x));
632 : }
633 49 : RgX_check_QX(x,"conjvec");
634 49 : av = avma;
635 49 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
636 49 : r = cleanroots(T,prec);
637 49 : z = cgetg(lx-2,t_COL);
638 182 : for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
639 49 : return gc_upto(av, z);
640 : }
641 :
642 0 : default:
643 0 : pari_err_TYPE("conjvec",x);
644 : return NULL; /* LCOV_EXCL_LINE */
645 : }
646 0 : return z;
647 : }
648 :
649 : /********************************************************************/
650 : /** **/
651 : /** ADDITION **/
652 : /** **/
653 : /********************************************************************/
654 : static GEN
655 24594801 : mkpadic_mod(GEN u, GEN p, GEN pd, long e, long d)
656 24594801 : { retmkpadic(modii(u, pd), icopy(p), icopy(pd), e, d); }
657 :
658 : /* x, y compatible PADIC, op = add or sub */
659 : static GEN
660 19974091 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
661 : {
662 19974091 : pari_sp av = avma;
663 : long d, e, r, rx, ry;
664 19974091 : GEN u, z, mod, pdx, pdy, ux, uy, p = padic_p(x);
665 : int swap;
666 :
667 19974091 : e = valp(x);
668 19974091 : r = valp(y); d = r-e;
669 19974091 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
670 19974091 : pdx = padic_pd(x); ux = padic_u(x);
671 19974091 : pdy = padic_pd(y); uy = padic_u(y);
672 19974091 : rx = precp(x);
673 19974091 : ry = precp(y);
674 19974091 : if (d) /* v(x) < v(y) */
675 : {
676 10704499 : r = d+ry; z = powiu(p,d);
677 10704623 : if (r < rx) mod = mulii(z,pdy); else { r = rx; mod = pdx; }
678 10704639 : z = mulii(z, uy);
679 10704595 : u = swap? op(z, ux): op(ux, z);
680 : }
681 : else
682 : {
683 : long c;
684 9269592 : if (ry < rx) { r=ry; mod = pdy; } else { r=rx; mod = pdx; }
685 9269592 : u = swap? op(uy, ux): op(ux, uy);
686 9271172 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
687 : {
688 138902 : set_avma(av); return zeropadic(p, e+r);
689 : }
690 9132081 : if (c)
691 : {
692 3351659 : mod = diviiexact(mod, powiu(p,c));
693 3351659 : r -= c;
694 3351659 : e += c;
695 : }
696 : }
697 19836644 : return gc_upto(av, mkpadic_mod(u, p, mod, e, r));
698 : }
699 : /* Rg_to_Fp(t_FRAC) without GC */
700 : static GEN
701 226838 : Q_to_Fp(GEN x, GEN p)
702 226838 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
703 : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
704 : static GEN
705 4763037 : addQp(GEN x, GEN y)
706 : {
707 4763037 : pari_sp av = avma;
708 4763037 : long d, r, e, vy = valp(y), py = precp(y);
709 4763037 : GEN q, mod, u, p = padic_p(y);
710 :
711 4763037 : e = Q_pvalrem(x, p, &x);
712 4763023 : d = vy - e; r = d + py;
713 4763023 : if (r <= 0) { set_avma(av); return gcopy(y); }
714 4761132 : mod = padic_pd(y);
715 4761132 : u = padic_u(y);
716 4761132 : if (d > 0)
717 : {
718 1376578 : q = powiu(p,d);
719 1376582 : mod = mulii(mod, q);
720 1376581 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
721 1376581 : u = addii(x, mulii(u, q));
722 : }
723 3384554 : else if (d < 0)
724 : {
725 405332 : q = powiu(p,-d);
726 405332 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
727 405332 : u = addii(u, mulii(x, q));
728 405332 : r = py; e = vy;
729 : }
730 : else
731 : {
732 : long c;
733 2979222 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
734 2979222 : u = addii(u, x);
735 2979243 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
736 : {
737 1204 : set_avma(av); return zeropadic(p,e+r);
738 : }
739 2978023 : if (c)
740 : {
741 970687 : mod = diviiexact(mod, powiu(p,c));
742 970686 : r -= c;
743 970686 : e += c;
744 : }
745 : }
746 4759940 : return gc_upto(av, mkpadic_mod(u, p, mod, e, r));
747 : }
748 :
749 : /* Mod(x,X) + Mod(y,X) */
750 : #define addsub_polmod_same addsub_polmod_scal
751 : /* Mod(x,X) +/- Mod(y,Y) */
752 : static GEN
753 7203 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
754 : {
755 7203 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
756 7203 : GEN z = cgetg(3,t_POLMOD);
757 7203 : long vx = varn(X), vy = varn(Y);
758 7203 : if (vx==vy) {
759 : pari_sp av;
760 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
761 14 : warn_coercion(X,Y,gel(z,1));
762 14 : gel(z,2) = gc_upto(av, gmod(op(x, y), gel(z,1))); return z;
763 : }
764 7189 : if (varncmp(vx, vy) < 0)
765 7189 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
766 : else
767 0 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
768 7189 : gel(z,2) = op(x, y); return z;
769 : }
770 : /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
771 : static GEN
772 13604332 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
773 13604332 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
774 :
775 : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
776 : static GEN
777 425209 : add_ser_scal(GEN y, GEN x)
778 : {
779 : long i, v, ly, vy;
780 : GEN z;
781 :
782 425209 : if (isrationalzero(x)) return gcopy(y);
783 398448 : ly = lg(y);
784 398448 : v = valser(y);
785 398448 : if (v < 3-ly) return gcopy(y);
786 : /* v + ly >= 3 */
787 398168 : if (v < 0)
788 : {
789 1176 : z = cgetg(ly,t_SER); z[1] = y[1];
790 3304 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
791 1176 : gel(z,i) = gadd(x,gel(y,i)); i++;
792 3087 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
793 1176 : return normalizeser(z);
794 : }
795 396992 : vy = varn(y);
796 396992 : if (v > 0)
797 : {
798 19670 : if (ser_isexactzero(y))
799 9506 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
800 10164 : y -= v; ly += v;
801 10164 : z = cgetg(ly,t_SER);
802 10164 : x = gcopy(x);
803 20762 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
804 : }
805 377322 : else if (ser_isexactzero(y)) return gcopy(y);
806 : else
807 : { /* v = 0, ly >= 3 */
808 377315 : z = cgetg(ly,t_SER);
809 377315 : x = gadd(x, gel(y,2));
810 377315 : i = 3;
811 : }
812 1648413 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
813 387479 : gel(z,2) = x;
814 387479 : z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
815 387479 : return gequal0(x)? normalizeser(z): z;
816 : }
817 : static long
818 8514516 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
819 : /* x,y t_SER in the same variable: x+y */
820 : static GEN
821 4257657 : ser_add(GEN x, GEN y)
822 : {
823 4257657 : long i, lx,ly, n = valser(y) - valser(x);
824 : GEN z;
825 4257657 : if (n < 0) { n = -n; swap(x,y); }
826 : /* valser(x) <= valser(y) */
827 4257657 : lx = _serprec(x);
828 4257657 : if (lx == 2) /* don't lose type information */
829 : {
830 798 : z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
831 798 : setvalser(z, valser(x)); return z;
832 : }
833 4256859 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
834 4256859 : if (n)
835 : {
836 107873 : if (n+2 > lx) return gcopy(x);
837 107173 : z = cgetg(ly,t_SER);
838 819459 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
839 506806 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
840 : } else {
841 4148986 : z = cgetg(ly,t_SER);
842 20825567 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
843 : }
844 4256159 : z[1] = x[1]; return normalizeser(z);
845 : }
846 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
847 : static GEN
848 8812479 : add_rfrac_scal(GEN y, GEN x)
849 : {
850 : pari_sp av;
851 : GEN n;
852 :
853 8812479 : if (isintzero(x)) return gcopy(y); /* frequent special case */
854 5138686 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
855 5138686 : return gc_upto(av, gred_rfrac_simple(n, gel(y,2)));
856 : }
857 :
858 : /* x "scalar", ty != t_MAT and nonscalar */
859 : static GEN
860 43070212 : add_scal(GEN y, GEN x, long ty)
861 : {
862 43070212 : switch(ty)
863 : {
864 38119042 : case t_POL: return RgX_Rg_add(y, x);
865 425181 : case t_SER: return add_ser_scal(y, x);
866 4480531 : case t_RFRAC: return add_rfrac_scal(y, x);
867 0 : case t_COL: return RgC_Rg_add(y, x);
868 45446 : case t_VEC:
869 45446 : if (isintzero(x)) return gcopy(y);
870 175 : break;
871 : }
872 187 : pari_err_TYPE2("+",x,y);
873 : return NULL; /* LCOV_EXCL_LINE */
874 : }
875 :
876 : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
877 : static GEN
878 15555600 : setfrac(GEN z, GEN a, GEN b)
879 : {
880 15555600 : gel(z,1) = icopy_avma(a, (pari_sp)z);
881 15555592 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
882 15555679 : set_avma((pari_sp)gel(z,2)); return z;
883 : }
884 : /* z <- a / (b*Q), (Q,a) = 1 */
885 : static GEN
886 14584195 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
887 : {
888 14584195 : GEN q = Qdivii(a, b); /* != 0 */
889 14584268 : if (typ(q) == t_INT)
890 : {
891 1925596 : gel(z,1) = gc_INT((pari_sp)Q, q);
892 1925596 : gel(z,2) = Q; return z;
893 : }
894 12658672 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
895 : }
896 : static GEN
897 41416191 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
898 : {
899 41416191 : GEN x1 = gel(x,1), x2 = gel(x,2);
900 41416191 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
901 41416191 : int s = cmpii(x2, y2);
902 :
903 : /* common denominator: (x1 op y1) / x2 */
904 41416184 : if (!s)
905 : {
906 15476388 : pari_sp av = avma;
907 15476388 : return gc_upto(av, Qdivii(op(x1, y1), x2));
908 : }
909 25939796 : z = cgetg(3, t_FRAC);
910 25939804 : if (s < 0)
911 : {
912 18085786 : Q = dvmdii(y2, x2, &r);
913 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
914 18085743 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
915 9849730 : delta = gcdii(x2,r);
916 : }
917 : else
918 : {
919 7854018 : Q = dvmdii(x2, y2, &r);
920 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
921 7854028 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
922 1505783 : delta = gcdii(y2,r);
923 : }
924 : /* delta = gcd(x2,y2) */
925 11355547 : if (equali1(delta))
926 : { /* numerator is nonzero */
927 8458534 : gel(z,1) = gc_INT((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
928 8458534 : gel(z,2) = mulii(x2,y2); return z;
929 : }
930 2897012 : x2 = diviiexact(x2,delta);
931 2897012 : y2 = diviiexact(y2,delta);
932 2897012 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
933 2897012 : q = dvmdii(n, delta, &r);
934 2897012 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
935 2642819 : r = gcdii(delta, r);
936 2642821 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
937 2642821 : return setfrac(z, n, mulii(mulii(x2, y2), delta));
938 : }
939 :
940 : /* assume x2, y2 are t_POLs in the same variable */
941 : static GEN
942 3017298 : add_rfrac(GEN x, GEN y)
943 : {
944 3017298 : pari_sp av = avma;
945 3017298 : GEN x1 = gel(x,1), x2 = gel(x,2);
946 3017298 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
947 :
948 3017298 : delta = RgX_gcd(x2,y2);
949 3017298 : if (!degpol(delta))
950 : {
951 658 : n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
952 658 : d = RgX_mul(x2, y2);
953 658 : return gc_upto(av, gred_rfrac_simple(n, d));
954 : }
955 3016640 : x2 = RgX_div(x2,delta);
956 3016640 : y2 = RgX_div(y2,delta);
957 3016640 : n = gadd(gmul(x1,y2), gmul(y1,x2));
958 3016640 : if (!signe(n))
959 : {
960 721446 : n = simplify_shallow(n);
961 721446 : if (isexactzero(n))
962 : {
963 721439 : if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
964 35 : return gc_upto(av, scalarpol(n, varn(x2)));
965 : }
966 7 : return gc_GEN(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
967 : }
968 2295194 : if (degpol(n) == 0)
969 1150589 : return gc_upto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
970 1144605 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
971 1144605 : if (isexactzero(r))
972 : {
973 : GEN z;
974 228024 : d = RgX_mul(x2, y2);
975 : /* "constant" denominator ? */
976 228024 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
977 228024 : return gc_upto(av, z);
978 : }
979 916581 : r = RgX_gcd(delta, r);
980 916581 : if (degpol(r))
981 : {
982 160593 : n = RgX_div(n, r);
983 160593 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
984 : }
985 : else
986 755988 : d = RgX_mul(gel(x,2), y2);
987 916581 : return gc_upto(av, gred_rfrac_simple(n, d));
988 : }
989 :
990 : GEN
991 6314982618 : gadd(GEN x, GEN y)
992 : {
993 6314982618 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
994 : pari_sp av;
995 : GEN z, p1;
996 :
997 6314982618 : if (tx == ty) switch(tx) /* shortcut to generic case */
998 : {
999 3088141457 : case t_INT: return addii(x,y);
1000 1751268570 : case t_REAL: return addrr(x,y);
1001 1244644 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1002 1244644 : z = cgetg(3,t_INTMOD);
1003 1244644 : if (X==Y || equalii(X,Y))
1004 1244630 : return add_intmod_same(z, X, gel(x,2), gel(y,2));
1005 14 : gel(z,1) = gcdii(X,Y);
1006 14 : warn_coercion(X,Y,gel(z,1));
1007 14 : av = avma; p1 = addii(gel(x,2),gel(y,2));
1008 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
1009 : }
1010 35236362 : case t_FRAC: return addsub_frac(x,y,addii);
1011 370367809 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1012 369950464 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1013 370717433 : if (isintzero(gel(z,2)))
1014 : {
1015 522851 : set_avma((pari_sp)(z+3));
1016 522851 : return gadd(gel(x,1),gel(y,1));
1017 : }
1018 369995789 : gel(z,1) = gadd(gel(x,1),gel(y,1));
1019 370169113 : return z;
1020 14531925 : case t_PADIC:
1021 14531925 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1022 14531620 : return addsub_pp(x,y, addii);
1023 672 : case t_QUAD: z = cgetg(4,t_QUAD);
1024 672 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1025 672 : gel(z,1) = ZX_copy(gel(x,1));
1026 672 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1027 672 : gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
1028 10213797 : case t_POLMOD:
1029 10213797 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1030 10206629 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
1031 7168 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
1032 33688964 : case t_FFELT: return FF_add(x,y);
1033 139510075 : case t_POL:
1034 139510075 : vx = varn(x);
1035 139510075 : vy = varn(y);
1036 139510075 : if (vx != vy) {
1037 4866601 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1038 1802362 : else return RgX_Rg_add(y, x);
1039 : }
1040 134643474 : return RgX_add(x, y);
1041 4252260 : case t_SER:
1042 4252260 : vx = varn(x);
1043 4252260 : vy = varn(y);
1044 4252260 : if (vx != vy) {
1045 28 : if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1046 21 : else return add_ser_scal(y, x);
1047 : }
1048 4252232 : return ser_add(x, y);
1049 4309199 : case t_RFRAC:
1050 4309199 : vx = varn(gel(x,2));
1051 4309199 : vy = varn(gel(y,2));
1052 4309199 : if (vx != vy) {
1053 1291901 : if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1054 538397 : else return add_rfrac_scal(y, x);
1055 : }
1056 3017298 : return add_rfrac(x,y);
1057 2118158 : case t_VEC:
1058 2118158 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1059 2118158 : return RgV_add(x,y);
1060 1137242 : case t_COL:
1061 1137242 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1062 1137242 : return RgC_add(x,y);
1063 671412 : case t_MAT:
1064 671412 : lx = lg(x);
1065 671412 : if (lg(y) != lx) pari_err_OP("+",x,y);
1066 671412 : if (lx == 1) return cgetg(1, t_MAT);
1067 671412 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1068 671405 : return RgM_add(x,y);
1069 :
1070 0 : default: pari_err_TYPE2("+",x,y);
1071 : }
1072 : /* tx != ty */
1073 862851135 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1074 :
1075 862851135 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1076 : {
1077 744576793 : case t_INT:
1078 744576793 : switch(ty)
1079 : {
1080 451349093 : case t_REAL: return addir(x,y);
1081 2141909 : case t_INTMOD:
1082 2141909 : z = cgetg(3, t_INTMOD);
1083 2141909 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1084 29062169 : case t_FRAC: z = cgetg(3,t_FRAC);
1085 29062164 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1086 29062138 : gel(z,2) = icopy(gel(y,2)); return z;
1087 255110668 : case t_COMPLEX: return addRc(x, y);
1088 5597067 : case t_PADIC:
1089 5597067 : if (!signe(x)) return gcopy(y);
1090 4533265 : return addQp(x,y);
1091 1421 : case t_QUAD: return addRq(x, y);
1092 1365414 : case t_FFELT: return FF_Z_add(y,x);
1093 : }
1094 :
1095 : case t_REAL:
1096 62181719 : switch(ty)
1097 : {
1098 14697610 : case t_FRAC:
1099 14697610 : if (!signe(gel(y,1))) return rcopy(x);
1100 14697610 : if (!signe(x))
1101 : {
1102 13590 : lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1103 13590 : return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1104 : }
1105 14684020 : av = avma; z = addir(gel(y,1), mulir(gel(y,2),x));
1106 14682315 : return gc_leaf(av, divri(z,gel(y,2)));
1107 47484039 : case t_COMPLEX: return addRc(x, y);
1108 70 : case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, realprec(x));
1109 :
1110 0 : default: pari_err_TYPE2("+",x,y);
1111 : }
1112 :
1113 17647 : case t_INTMOD:
1114 17647 : switch(ty)
1115 : {
1116 17507 : case t_FRAC: { GEN X = gel(x,1);
1117 17507 : z = cgetg(3, t_INTMOD);
1118 17507 : p1 = Fp_div(gel(y,1), gel(y,2), X);
1119 17507 : return add_intmod_same(z, X, p1, gel(x,2));
1120 : }
1121 14 : case t_FFELT:
1122 14 : if (!equalii(gel(x,1),FF_p_i(y)))
1123 0 : pari_err_OP("+",x,y);
1124 14 : return FF_Z_add(y,gel(x,2));
1125 91 : case t_COMPLEX: return addRc(x, y);
1126 0 : case t_PADIC: { GEN X = gel(x,1);
1127 0 : z = cgetg(3, t_INTMOD);
1128 0 : return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1129 : }
1130 35 : case t_QUAD: return addRq(x, y);
1131 : }
1132 :
1133 : case t_FRAC:
1134 10686204 : switch (ty)
1135 : {
1136 10454962 : case t_COMPLEX: return addRc(x, y);
1137 229771 : case t_PADIC:
1138 229771 : if (!signe(gel(x,1))) return gcopy(y);
1139 229771 : return addQp(x,y);
1140 133 : case t_QUAD: return addRq(x, y);
1141 1337 : case t_FFELT: return FF_Q_add(y, x);
1142 : }
1143 :
1144 : case t_FFELT:
1145 1 : pari_err_TYPE2("+",x,y);
1146 :
1147 35 : case t_COMPLEX:
1148 35 : switch(ty)
1149 : {
1150 28 : case t_PADIC:
1151 28 : return Zp_nosquare_m1(padic_p(y))? addRc(y, x): addTp(x, y);
1152 7 : case t_QUAD:
1153 7 : lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1154 7 : return gequal0(y)? gcopy(x): addqf(y, x, lx);
1155 : }
1156 :
1157 : case t_PADIC: /* ty == t_QUAD */
1158 6 : return (kro_quad(y, padic_p(x)) == -1)? addRq(x, y): addTp(y, x);
1159 : }
1160 : /* tx < ty, !is_const_t(y) */
1161 48758568 : switch(ty)
1162 : {
1163 7948 : case t_MAT:
1164 7948 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1165 7948 : if (isrationalzero(x)) return gcopy(y);
1166 7871 : return RgM_Rg_add(y, x);
1167 186499 : case t_COL:
1168 186499 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1169 186499 : return RgC_Rg_add(y, x);
1170 2417220 : case t_POLMOD: /* is_const_t(tx) in this case */
1171 2417220 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1172 : }
1173 46146901 : if (is_scalar_t(tx)) {
1174 43078125 : if (tx == t_POLMOD)
1175 : {
1176 118613 : vx = varn(gel(x,1));
1177 118613 : vy = gvar(y);
1178 118613 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1179 : else
1180 88485 : if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1181 30667 : return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1182 : }
1183 42959512 : return add_scal(y, x, ty);
1184 : }
1185 : /* x and y are not scalars, ty != t_MAT */
1186 3068766 : vx = gvar(x);
1187 3068777 : vy = gvar(y);
1188 3068777 : if (vx != vy) { /* x or y is treated as a scalar */
1189 22773 : if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1190 32445 : return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1191 32445 : : add_scal(y, x, ty);
1192 : }
1193 : /* vx = vy */
1194 3046004 : switch(tx)
1195 : {
1196 3045500 : case t_POL:
1197 3045500 : switch (ty)
1198 : {
1199 5453 : case t_SER:
1200 5453 : if (lg(x) == 2) return gcopy(y);
1201 5432 : i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1202 5432 : i = lg(y) + valser(y) - i;
1203 5432 : if (i < 3) return gcopy(y);
1204 5425 : p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1205 5425 : settyp(p1, t_VECSMALL); /* p1 left on stack */
1206 5425 : return y;
1207 :
1208 3040047 : case t_RFRAC: return add_rfrac_scal(y, x);
1209 : }
1210 0 : break;
1211 :
1212 504 : case t_SER:
1213 504 : if (ty == t_RFRAC)
1214 : {
1215 : long vn, vd;
1216 504 : av = avma;
1217 504 : vn = gval(gel(y,1), vy);
1218 504 : vd = RgX_valrem_inexact(gel(y,2), NULL);
1219 504 : if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
1220 :
1221 504 : l = lg(x) + valser(x) - (vn - vd);
1222 504 : if (l < 3) { set_avma(av); return gcopy(x); }
1223 504 : return gc_upto(av, gadd(x, rfrac_to_ser_i(y, l)));
1224 : }
1225 0 : break;
1226 : }
1227 0 : pari_err_TYPE2("+",x,y);
1228 : return NULL; /* LCOV_EXCL_LINE */
1229 : }
1230 :
1231 : GEN
1232 270841641 : gaddsg(long x, GEN y)
1233 : {
1234 270841641 : long ty = typ(y);
1235 : GEN z;
1236 :
1237 270841641 : switch(ty)
1238 : {
1239 124845856 : case t_INT: return addsi(x,y);
1240 119317089 : case t_REAL: return addsr(x,y);
1241 12047 : case t_INTMOD:
1242 12047 : z = cgetg(3, t_INTMOD);
1243 12047 : return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1244 15251863 : case t_FRAC: z = cgetg(3,t_FRAC);
1245 15251863 : gel(z,1) = gc_INT((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1246 15251863 : gel(z,2) = icopy(gel(y,2)); return z;
1247 8270647 : case t_COMPLEX:
1248 8270647 : z = cgetg(3, t_COMPLEX);
1249 8270647 : gel(z,1) = gaddsg(x, gel(y,1));
1250 8270647 : gel(z,2) = gcopy(gel(y,2)); return z;
1251 :
1252 3144139 : default: return gadd(stoi(x), y);
1253 : }
1254 : }
1255 :
1256 : GEN
1257 3271588 : gsubsg(long x, GEN y)
1258 : {
1259 : GEN z, a, b;
1260 : pari_sp av;
1261 :
1262 3271588 : switch(typ(y))
1263 : {
1264 277975 : case t_INT: return subsi(x,y);
1265 1247433 : case t_REAL: return subsr(x,y);
1266 56 : case t_INTMOD:
1267 56 : z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1268 56 : return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1269 732541 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1270 732541 : gel(z,1) = gc_INT((pari_sp)z, subii(mulis(b,x), a));
1271 732540 : gel(z,2) = icopy(gel(y,2)); return z;
1272 974607 : case t_COMPLEX:
1273 974607 : z = cgetg(3, t_COMPLEX);
1274 974607 : gel(z,1) = gsubsg(x, gel(y,1));
1275 974607 : gel(z,2) = gneg(gel(y,2)); return z;
1276 : }
1277 38976 : av = avma;
1278 38976 : return gc_upto(av, gadd(stoi(x), gneg_i(y)));
1279 : }
1280 :
1281 : /********************************************************************/
1282 : /** **/
1283 : /** SUBTRACTION **/
1284 : /** **/
1285 : /********************************************************************/
1286 :
1287 : GEN
1288 2933498489 : gsub(GEN x, GEN y)
1289 : {
1290 2933498489 : long tx = typ(x), ty = typ(y);
1291 : pari_sp av;
1292 : GEN z;
1293 2933498489 : if (tx == ty) switch(tx) /* shortcut to generic case */
1294 : {
1295 2093679992 : case t_INT: return subii(x,y);
1296 619064530 : case t_REAL: return subrr(x,y);
1297 1031415 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1298 1031415 : z = cgetg(3,t_INTMOD);
1299 1031415 : if (X==Y || equalii(X,Y))
1300 1031401 : return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1301 14 : gel(z,1) = gcdii(X,Y);
1302 14 : warn_coercion(X,Y,gel(z,1));
1303 14 : av = avma; p1 = subii(gel(x,2),gel(y,2));
1304 14 : gel(z,2) = gc_INT(av, modii(p1, gel(z,1))); return z;
1305 : }
1306 6179849 : case t_FRAC: return addsub_frac(x,y, subii);
1307 104653897 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1308 104687892 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1309 104602035 : if (isintzero(gel(z,2)))
1310 : {
1311 21322 : set_avma((pari_sp)(z+3));
1312 21322 : return gsub(gel(x,1),gel(y,1));
1313 : }
1314 104580636 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1315 104570242 : return z;
1316 5442581 : case t_PADIC:
1317 5442581 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("+",x,y);
1318 5442581 : return addsub_pp(x,y, subii);
1319 168 : case t_QUAD: z = cgetg(4,t_QUAD);
1320 168 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1321 168 : gel(z,1) = ZX_copy(gel(x,1));
1322 168 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1323 168 : gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1324 949851 : case t_POLMOD:
1325 949851 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1326 949816 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1327 35 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1328 7934937 : case t_FFELT: return FF_sub(x,y);
1329 20326495 : case t_POL: {
1330 20326495 : long vx = varn(x);
1331 20326495 : long vy = varn(y);
1332 20326495 : if (vx != vy) {
1333 1154691 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1334 86024 : else return Rg_RgX_sub(x, y);
1335 : }
1336 19171804 : return RgX_sub(x, y);
1337 : }
1338 324365 : case t_VEC:
1339 324365 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1340 324365 : return RgV_sub(x,y);
1341 3605021 : case t_COL:
1342 3605021 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1343 3605021 : return RgC_sub(x,y);
1344 254127 : case t_MAT: {
1345 254127 : long lx = lg(x);
1346 254127 : if (lg(y) != lx) pari_err_OP("+",x,y);
1347 254129 : if (lx == 1) return cgetg(1, t_MAT);
1348 173911 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1349 173914 : return RgM_sub(x,y);
1350 : }
1351 2583358 : case t_RFRAC: case t_SER: break;
1352 :
1353 0 : default: pari_err_TYPE2("+",x,y);
1354 : }
1355 71668777 : av = avma;
1356 71668777 : return gc_upto(av, gadd(x,gneg_i(y)));
1357 : }
1358 :
1359 : /********************************************************************/
1360 : /** **/
1361 : /** MULTIPLICATION **/
1362 : /** **/
1363 : /********************************************************************/
1364 : static GEN
1365 326494 : mul_ser_scal(GEN x, GEN t)
1366 : {
1367 326494 : if (isexactzero(t)) return gmul(Rg_get_0(x), t);
1368 323246 : if (isint1(t)) return gcopy(x);
1369 268695 : if (ser_isexactzero(x))
1370 : {
1371 378 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gmul(gel(x,2), t), varn(x), 1);
1372 378 : setvalser(z, valser(x)); return z;
1373 : }
1374 3465952 : pari_APPLY_ser(gmul(gel(x,i), t));
1375 : }
1376 : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1377 : * [n/d a valid RFRAC] */
1378 : static GEN
1379 10454352 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1380 : {
1381 10454352 : pari_sp av = avma;
1382 : GEN z;
1383 :
1384 10454352 : switch(typ(x))
1385 : {
1386 21 : case t_PADIC:
1387 21 : n = gmul(n, x);
1388 21 : d = gcvtop(d, padic_p(x), signe(padic_u(x))? precp(x): 1);
1389 21 : return gc_upto(av, gdiv(n,d));
1390 :
1391 952 : case t_INTMOD: case t_POLMOD:
1392 952 : n = gmul(n, x);
1393 952 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1394 952 : return gc_upto(av, gdiv(n,d));
1395 : }
1396 10453379 : z = gred_rfrac(x, d);
1397 10453379 : n = simplify_shallow(n);
1398 10453379 : if (typ(z) == t_RFRAC)
1399 : {
1400 7928812 : n = gmul(gel(z,1), n);
1401 7928812 : d = gel(z,2);
1402 7928812 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1403 0 : z = RgX_Rg_div(n, d);
1404 : else
1405 7928812 : z = gred_rfrac_simple(n, d);
1406 : }
1407 : else
1408 2524567 : z = gmul(z, n);
1409 10453379 : return gc_upto(av, z);
1410 : }
1411 : static GEN
1412 155040995 : mul_scal(GEN y, GEN x, long ty)
1413 : {
1414 155040995 : switch(ty)
1415 : {
1416 146307301 : case t_POL:
1417 146307301 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1418 134234442 : return RgX_Rg_mul(y, x);
1419 318402 : case t_SER: return mul_ser_scal(y, x);
1420 8415269 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1421 14 : case t_QFB:
1422 14 : if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1423 : }
1424 12 : pari_err_TYPE2("*",x,y);
1425 : return NULL; /* LCOV_EXCL_LINE */
1426 : }
1427 : static GEN
1428 8027393 : mul_self_scal(GEN x, GEN y)
1429 642931188 : { pari_APPLY_same(gmul(y,gel(x,i))); }
1430 :
1431 : static GEN
1432 161432 : mul_gen_rfrac(GEN X, GEN Y)
1433 : {
1434 161432 : GEN y1 = gel(Y,1), y2 = gel(Y,2);
1435 161432 : long vx = gvar(X), vy = varn(y2);
1436 167193 : return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1437 5761 : gred_rfrac_simple(gmul(y1,X), y2);
1438 : }
1439 : /* (x1/x2) * (y1/y2) */
1440 : static GEN
1441 7910397 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1442 : {
1443 : GEN z, X, Y;
1444 7910397 : pari_sp av = avma;
1445 :
1446 7910397 : X = gred_rfrac(x1, y2);
1447 7910397 : Y = gred_rfrac(y1, x2);
1448 7910397 : if (typ(X) == t_RFRAC)
1449 : {
1450 6631301 : if (typ(Y) == t_RFRAC) {
1451 6565144 : x1 = gel(X,1);
1452 6565144 : x2 = gel(X,2);
1453 6565144 : y1 = gel(Y,1);
1454 6565144 : y2 = gel(Y,2);
1455 6565144 : z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1456 : } else
1457 66157 : z = mul_gen_rfrac(Y, X);
1458 : }
1459 1279096 : else if (typ(Y) == t_RFRAC)
1460 95275 : z = mul_gen_rfrac(X, Y);
1461 : else
1462 1183821 : z = gmul(X, Y);
1463 7910397 : return gc_upto(av, z);
1464 : }
1465 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1466 : static GEN
1467 262365 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1468 : {
1469 262365 : pari_sp av = avma;
1470 262365 : GEN X = gred_rfrac(x1, y2);
1471 262365 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1472 : {
1473 255785 : x2 = RgX_mul(gel(X,2), x2);
1474 255785 : x1 = gel(X,1);
1475 : }
1476 : else
1477 6580 : x1 = X;
1478 262365 : return gc_upto(av, gred_rfrac_simple(x1, x2));
1479 : }
1480 :
1481 : /* Mod(y, Y) * x, assuming x scalar */
1482 : static GEN
1483 3288167 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1484 3288167 : { retmkpolmod(isrationalzero(x)? gen_0: gmul(x,y), RgX_copy(Y)); }
1485 :
1486 : /* cf mulqq */
1487 : static GEN
1488 5861735 : quad_polmod_mul(GEN T, GEN x, GEN y)
1489 : {
1490 5861735 : GEN b = gel(T,3), c = gel(T,2);
1491 5861735 : GEN ux = gel(x,2), vx = gel(x,3), uy = gel(y,2), vy = gel(y,3);
1492 : GEN z, s, U, V, E, F;
1493 : pari_sp av, av2;
1494 5861735 : z = cgetg(4, t_POL); z[1] = x[1]; av = avma;
1495 5861735 : U = gmul(ux, uy);
1496 5861735 : V = gmul(vx, vy); s = gmul(c, V);
1497 5861735 : E = gmul(gadd(ux, vx), gadd(uy, vy));
1498 5861735 : if (typ(b) != t_INT) F = gadd(U, gmul(gaddgs(b, 1), V)); /* = U + (b+1) V */
1499 : else
1500 : { /* minor optimization */
1501 5861714 : if (!signe(b)) F = gadd(U, V); /* b = 0 */
1502 4286487 : else if (is_pm1(b)) /* b = 1 or -1 */
1503 4285829 : F = signe(b) < 0? U: gadd(U, gmul2n(V,1));
1504 : else /* |b| > 1 */
1505 658 : F = gadd(U, gmul(addis(b, 1), V));
1506 : }
1507 5861735 : av2 = avma;
1508 5861735 : gel(z,2) = gsub(U, s);
1509 5861735 : gel(z,3) = gsub(E, F);
1510 5861735 : gc_slice_unsafe(av, av2, z+2, 2);
1511 5861735 : return normalizepol_lg(z,4);
1512 : }
1513 : /* Mod(x,T) * Mod(y,T) */
1514 : static GEN
1515 8646527 : mul_polmod_same(GEN T, GEN x, GEN y)
1516 : {
1517 8646527 : GEN z = cgetg(3,t_POLMOD), a;
1518 8646527 : long v = varn(T), lx = lg(x), ly = lg(y);
1519 8646527 : gel(z,1) = RgX_copy(T);
1520 : /* x * y mod T optimised */
1521 8646527 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1522 7905771 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1523 1662002 : a = gmul(x, y);
1524 : else
1525 : {
1526 6984525 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1527 5856387 : a = quad_polmod_mul(T, x, y);
1528 : else
1529 1128138 : a = RgXQ_mul(x, y, gel(z,1));
1530 : }
1531 8646527 : gel(z,2) = a; return z;
1532 : }
1533 : static GEN
1534 426147 : sqr_polmod(GEN T, GEN x)
1535 : {
1536 426147 : GEN a, z = cgetg(3,t_POLMOD);
1537 426147 : gel(z,1) = RgX_copy(T);
1538 426147 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1539 35407 : a = gsqr(x);
1540 : else
1541 : {
1542 390740 : pari_sp av = avma;
1543 390740 : a = RgXQ_sqr(x, gel(z,1));
1544 390740 : a = gc_upto(av, a);
1545 : }
1546 426147 : gel(z,2) = a; return z;
1547 : }
1548 : /* Mod(x,X) * Mod(y,Y) */
1549 : static GEN
1550 2674945 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1551 : {
1552 2674945 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1553 2674945 : long vx = varn(X), vy = varn(Y);
1554 2674945 : GEN z = cgetg(3,t_POLMOD);
1555 :
1556 2674945 : if (vx==vy) {
1557 : pari_sp av;
1558 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
1559 14 : warn_coercion(X,Y,gel(z,1));
1560 14 : gel(z,2) = gc_upto(av, gmod(gmul(x, y), gel(z,1)));
1561 14 : return z;
1562 : }
1563 2674931 : if (varncmp(vx, vy) < 0)
1564 410662 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1565 : else
1566 2264269 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1567 2674931 : gel(z,2) = gmul(x, y); return z;
1568 : }
1569 :
1570 : #if 0 /* used by 3M only */
1571 : /* set z = x+y and return 1 if x,y have the same sign
1572 : * set z = x-y and return 0 otherwise */
1573 : static int
1574 : did_add(GEN x, GEN y, GEN *z)
1575 : {
1576 : long tx = typ(x), ty = typ(y);
1577 : if (tx == ty) switch(tx)
1578 : {
1579 : case t_INT: *z = addii(x,y); return 1;
1580 : case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1581 : case t_REAL:
1582 : if (signe(x) == -signe(y))
1583 : { *z = subrr(x,y); return 0; }
1584 : else
1585 : { *z = addrr(x,y); return 1; }
1586 : }
1587 : if (tx == t_REAL) switch(ty)
1588 : {
1589 : case t_INT:
1590 : if (signe(x) == -signe(y))
1591 : { *z = subri(x,y); return 0; }
1592 : else
1593 : { *z = addri(x,y); return 1; }
1594 : case t_FRAC:
1595 : if (signe(x) == -signe(gel(y,1)))
1596 : { *z = gsub(x,y); return 0; }
1597 : else
1598 : { *z = gadd(x,y); return 1; }
1599 : }
1600 : else if (ty == t_REAL) switch(tx)
1601 : {
1602 : case t_INT:
1603 : if (signe(x) == -signe(y))
1604 : { *z = subir(x,y); return 0; }
1605 : else
1606 : { *z = addir(x,y); return 1; }
1607 : case t_FRAC:
1608 : if (signe(gel(x,1)) == -signe(y))
1609 : { *z = gsub(x,y); return 0; }
1610 : else
1611 : { *z = gadd(x,y); return 1; }
1612 : }
1613 : *z = gadd(x,y); return 1;
1614 : }
1615 : #endif
1616 : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1617 : static GEN
1618 11900098 : mulcIR(GEN x, GEN y)
1619 : {
1620 11900098 : GEN z = cgetg(3,t_COMPLEX);
1621 11900238 : pari_sp av = avma;
1622 11900238 : gel(z,1) = gc_upto(av, gneg(gmul(y,gel(x,2))));
1623 11900378 : gel(z,2) = gmul(y, gel(x,1));
1624 11899800 : return z;
1625 :
1626 : }
1627 : /* x,y COMPLEX */
1628 : static GEN
1629 286251189 : mulcc(GEN x, GEN y)
1630 : {
1631 286251189 : GEN xr = gel(x,1), xi = gel(x,2);
1632 286251189 : GEN yr = gel(y,1), yi = gel(y,2);
1633 : GEN p1, p2, p3, p4, z;
1634 : pari_sp tetpil, av;
1635 :
1636 286251189 : if (isintzero(xr))
1637 : {
1638 15653432 : if (isintzero(yr)) {
1639 7122788 : av = avma;
1640 7122788 : return gc_upto(av, gneg(gmul(xi,yi)));
1641 : }
1642 8530510 : return mulcIR(y, xi);
1643 : }
1644 270587660 : if (isintzero(yr)) return mulcIR(x, yi);
1645 :
1646 267203273 : z = cgetg(3,t_COMPLEX); av = avma;
1647 : #if 0
1648 : /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1649 : * e.g. xr + xi if exponents differ */
1650 : if (did_add(xr, xi, &p3))
1651 : {
1652 : if (did_add(yr, yi, &p4)) {
1653 : /* R = xr*yr - xi*yi
1654 : * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1655 : p1 = gmul(xr,yr);
1656 : p2 = gmul(xi,yi); p2 = gneg(p2);
1657 : p3 = gmul(p3, p4);
1658 : p4 = gsub(p2, p1);
1659 : } else {
1660 : /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1661 : * I = xr*yi + xi*yr */
1662 : p1 = gmul(p3,p4);
1663 : p3 = gmul(xr,yi);
1664 : p4 = gmul(xi,yr);
1665 : p2 = gsub(p3, p4);
1666 : }
1667 : } else {
1668 : if (did_add(yr, yi, &p4)) {
1669 : /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1670 : * I = xr*yi +xi*yr */
1671 : p1 = gmul(p3,p4);
1672 : p3 = gmul(xr,yi);
1673 : p4 = gmul(xi,yr);
1674 : p2 = gsub(p4, p3);
1675 : } else {
1676 : /* R = xr*yr - xi*yi
1677 : * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1678 : p3 = gneg( gmul(p3, p4) );
1679 : p1 = gmul(xr,yr);
1680 : p2 = gmul(xi,yi);
1681 : p4 = gadd(p1, p2);
1682 :
1683 : p2 = gneg(p2);
1684 : }
1685 : }
1686 : tetpil = avma;
1687 : gel(z,1) = gadd(p1,p2);
1688 : gel(z,2) = gadd(p3,p4);
1689 : #else
1690 267205534 : if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1691 : { /* 3M formula */
1692 558740 : p3 = addii(xr,xi);
1693 558740 : p4 = addii(yr,yi);
1694 558740 : p1 = mulii(xr,yr);
1695 558740 : p2 = mulii(xi,yi);
1696 558740 : p3 = mulii(p3,p4);
1697 558740 : p4 = addii(p2,p1);
1698 558740 : tetpil = avma;
1699 558740 : gel(z,1) = subii(p1,p2);
1700 558740 : gel(z,2) = subii(p3,p4);
1701 558740 : if (!signe(gel(z,2)))
1702 113260 : return gc_INT((pari_sp)(z+3), gel(z,1));
1703 : }
1704 : else
1705 : { /* naive 4M formula: avoid all loss of accuracy */
1706 266646794 : p1 = gmul(xr,yr);
1707 266633283 : p2 = gmul(xi,yi);
1708 266623414 : p3 = gmul(xr,yi);
1709 266623917 : p4 = gmul(xi,yr);
1710 266614278 : tetpil = avma;
1711 266614278 : gel(z,1) = gsub(p1,p2);
1712 266491737 : gel(z,2) = gadd(p3,p4);
1713 266498971 : if (isintzero(gel(z,2)))
1714 : {
1715 50701 : cgiv(gel(z,2));
1716 50701 : return gc_upto((pari_sp)(z+3), gel(z,1));
1717 : }
1718 : }
1719 : #endif
1720 :
1721 266889348 : gc_slice_unsafe(av,tetpil, z+1,2); return z;
1722 : }
1723 : /* x,y PADIC */
1724 : static GEN
1725 17505614 : mulpp(GEN x, GEN y) {
1726 17505614 : GEN pd, p = padic_p(x), ux = padic_u(x), uy = padic_u(y);
1727 17505614 : long e = valp(x) + valp(y), dx, dy;
1728 :
1729 17505614 : if (!equalii(p, padic_p(y))) pari_err_OP("*",x,y);
1730 17505537 : if (!signe(ux) || !signe(uy)) return zeropadic(p, e);
1731 16931598 : dx = precp(x); dy = precp(y);
1732 16931598 : if (dx > dy) { pd = padic_pd(y); dx = dy; } else pd = padic_pd(x);
1733 16931598 : retmkpadic(Fp_mul(ux, uy, pd), icopy(p), icopy(pd), e, dx);
1734 : }
1735 : /* x,y QUAD, w^2 = - b w - c, where b = 0 or -1
1736 : * (ux + vx w)(uy + vy w) = ux uy - c vx vy + (ux vy + uy vx - b vx vy) w */
1737 : static GEN
1738 1127 : mulqq(GEN x, GEN y)
1739 : {
1740 1127 : GEN T = gel(x,1), b = gel(T,3), c = gel(T,2);
1741 1127 : GEN ux = gel(x,2), vx = gel(x,3), uy = gel(y,2), vy = gel(y,3);
1742 : GEN z, s, U, V, E, F;
1743 : pari_sp av, av2;
1744 1127 : z = cgetg(4, t_QUAD), gel(z,1) = ZX_copy(T); av = avma;
1745 1127 : U = gmul(ux, uy);
1746 1127 : V = gmul(vx, vy); s = gmul(c, V);
1747 1127 : E = gmul(gadd(ux, vx), gadd(uy, vy));
1748 1127 : F = signe(b)? U: gadd(U, V);
1749 : /* E - F = ux vy + uy vx - b vx vy */
1750 1127 : av2 = avma;
1751 1127 : gel(z,2) = gsub(U, s);
1752 1127 : gel(z,3) = gsub(E, F);
1753 1127 : gc_slice_unsafe(av, av2, z+2, 2); return z;
1754 : }
1755 :
1756 : GEN
1757 15677916 : mulcxI(GEN x)
1758 : {
1759 : GEN z;
1760 15677916 : switch(typ(x))
1761 : {
1762 1009782 : case t_INT:
1763 1009782 : if (!signe(x)) return gen_0;
1764 : case t_REAL: case t_FRAC:
1765 2623177 : return mkcomplex(gen_0, x);
1766 13054452 : case t_COMPLEX:
1767 13054452 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1768 13041042 : z = cgetg(3,t_COMPLEX);
1769 13041313 : gel(z,1) = gneg(gel(x,2));
1770 13042017 : gel(z,2) = gel(x,1); return z;
1771 112 : default:
1772 112 : return gmul(gen_I(), x);
1773 : }
1774 : }
1775 : GEN
1776 3166425 : mulcxmI(GEN x)
1777 : {
1778 : GEN z;
1779 3166425 : switch(typ(x))
1780 : {
1781 546 : case t_INT:
1782 546 : if (!signe(x)) return gen_0;
1783 : case t_REAL: case t_FRAC:
1784 141235 : return mkcomplex(gen_0, gneg(x));
1785 2822216 : case t_COMPLEX:
1786 2822216 : if (isintzero(gel(x,1))) return gel(x,2);
1787 2727231 : z = cgetg(3,t_COMPLEX);
1788 2727198 : gel(z,1) = gel(x,2);
1789 2727198 : gel(z,2) = gneg(gel(x,1)); return z;
1790 202974 : default:
1791 202974 : return gmul(mkcomplex(gen_0, gen_m1), x);
1792 : }
1793 : }
1794 : /* x * I^k */
1795 : GEN
1796 5746459 : mulcxpowIs(GEN x, long k)
1797 : {
1798 5746459 : switch (k & 3)
1799 : {
1800 1369672 : case 1: return mulcxI(x);
1801 1362757 : case 2: return gneg(x);
1802 1583932 : case 3: return mulcxmI(x);
1803 : }
1804 1430098 : return x;
1805 : }
1806 :
1807 : /* Convert t_POL y to t_SER of length l and valuation e then normalize. Must
1808 : * not be called when y is an exact zero carrying type information, see comment
1809 : * in gmul(t_SER, t_SER): the valuation may incorrectly be altered by
1810 : * normalizeser in this case. Variant of RgX_to_ser_i with a different
1811 : * meaning for e. */
1812 : static GEN
1813 6631815 : pol2ser(GEN y, long l, long e)
1814 : {
1815 6631815 : long i, ly = lg(y);
1816 6631815 : GEN z = cgetg(l, t_SER);
1817 6631815 : z[1] = evalvalser(e) | evalvarn(varn(y)) | evalsigne(1);
1818 6631815 : if (ly >= l) {
1819 30196627 : for (i = 2; i < l; i++) gel(z,i) = gel(y,i);
1820 : } else {
1821 339909 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1822 240644 : for ( ; i < l; i++) gel(z,i) = gen_0;
1823 : }
1824 6631815 : return normalizeser(z);
1825 : }
1826 : /* assume typ(x) = t_VEC */
1827 : static int
1828 98 : is_ext_qfr(GEN x)
1829 91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
1830 189 : && typ(gel(x,2)) == t_REAL; }
1831 :
1832 : GEN
1833 8880109104 : gmul(GEN x, GEN y)
1834 : {
1835 : long tx, ty, lx, ly, vx, vy, i, l;
1836 : pari_sp av, tetpil;
1837 : GEN z, p1;
1838 :
1839 8880109104 : if (x == y) return gsqr(x);
1840 7918065441 : tx = typ(x); ty = typ(y);
1841 7918065441 : if (tx == ty) switch(tx)
1842 : {
1843 3556722357 : case t_INT: return mulii(x,y);
1844 2006349005 : case t_REAL: return mulrr(x,y);
1845 1451928 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1846 1451928 : z = cgetg(3,t_INTMOD);
1847 1451928 : if (X==Y || equalii(X,Y))
1848 1451886 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1849 42 : gel(z,1) = gcdii(X,Y);
1850 42 : warn_coercion(X,Y,gel(z,1));
1851 42 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1852 42 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
1853 : }
1854 23558271 : case t_FRAC:
1855 : {
1856 23558271 : GEN x1 = gel(x,1), x2 = gel(x,2);
1857 23558271 : GEN y1 = gel(y,1), y2 = gel(y,2);
1858 23558271 : z=cgetg(3,t_FRAC);
1859 23558270 : p1 = gcdii(x1, y2);
1860 23558270 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1861 23558270 : p1 = gcdii(x2, y1);
1862 23558269 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1863 23558269 : tetpil = avma;
1864 23558269 : gel(z,2) = mulii(x2,y2);
1865 23558267 : gel(z,1) = mulii(x1,y1);
1866 23558269 : fix_frac_if_int_GC(z,tetpil); return z;
1867 : }
1868 276043619 : case t_COMPLEX: return mulcc(x, y);
1869 12511146 : case t_PADIC: return mulpp(x, y);
1870 882 : case t_QUAD:
1871 882 : if (!ZX_equal(gel(x,1), gel(y,1))) pari_err_OP("*",x,y);
1872 882 : return mulqq(x, y);
1873 25760666 : case t_FFELT: return FF_mul(x, y);
1874 11068373 : case t_POLMOD:
1875 11068373 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1876 8393428 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1877 2674945 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1878 132158821 : case t_POL:
1879 132158821 : vx = varn(x);
1880 132158821 : vy = varn(y);
1881 132158821 : if (vx != vy) {
1882 46850710 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1883 27719338 : else return RgX_Rg_mul(y, x);
1884 : }
1885 85308111 : return RgX_mul(x, y);
1886 :
1887 5243415 : case t_SER:
1888 : {
1889 : long e;
1890 5243415 : vx = varn(x);
1891 5243415 : vy = varn(y);
1892 5243415 : if (vx != vy) {
1893 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1894 3675 : return mul_ser_scal(y, x);
1895 : }
1896 5239740 : e = valser(x) + valser(y);
1897 5239740 : lx = minss(lg(x), lg(y));
1898 5239740 : if (lx == 2) return zeroser(vx, e);
1899 : /* Must handle the case of exact zeros here. If x or y is t^e (z+O(t))
1900 : * where z is an exact zero giving type information, then its valser is
1901 : * v+1 and their product has the same shape but we mustn't call
1902 : * normalizeser() since the result is already normalized (see "dangerous
1903 : * case" there). It would add 1 to the valuation again.
1904 : *
1905 : * But x,y may be zero divisors, Mod(2,4) + O(t) say, and the
1906 : * product again has the above shape but this time we must normalize.
1907 : * There is no way to distinguish between the two situations in
1908 : * pol2ser() from the product of the polynomial parts only */
1909 5239726 : if (lx == 3)
1910 : {
1911 : long s;
1912 : GEN c;
1913 26509 : z = cgetg(3, t_SER); c = gmul(gel(x,2), gel(y,2));
1914 26509 : s = !gequal0(c);
1915 26509 : if (!s && signe(x) && signe(y) && isexactzero(c)) e++;
1916 26509 : z[1] = evalsigne(s) | evalvarn(vx) | evalvalser(e);
1917 26509 : gel(z,2) = c; return z;
1918 : }
1919 5213217 : av = avma;
1920 5213217 : x = ser2pol_i(x, lx);
1921 5213217 : y = ser2pol_i(y, lx);
1922 5213217 : y = RgXn_mul(x, y, lx-2);
1923 5213217 : return gc_GEN(av, pol2ser(y, lx, e));
1924 : }
1925 42 : case t_VEC:
1926 42 : if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
1927 : /* fall through, handle extended t_QFB */
1928 198 : case t_QFB: return qfbcomp(x,y);
1929 6723146 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1930 4129066 : case t_MAT: return RgM_mul(x, y);
1931 :
1932 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1933 1631 : z = cgetg_copy(x, &l);
1934 1631 : if (l != lg(y)) break;
1935 17325 : for (i=1; i<l; i++)
1936 : {
1937 15694 : long yi = y[i];
1938 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1939 15694 : z[i] = x[yi];
1940 : }
1941 1631 : return z;
1942 :
1943 0 : default:
1944 0 : pari_err_TYPE2("*",x,y);
1945 : }
1946 : /* tx != ty */
1947 1858165518 : if (is_const_t(ty) && is_const_t(tx)) {
1948 1694349241 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1949 1694349241 : switch(tx) {
1950 1562728560 : case t_INT:
1951 1562728560 : switch(ty)
1952 : {
1953 1015674038 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1954 1686310 : case t_INTMOD:
1955 1686310 : z = cgetg(3, t_INTMOD);
1956 1686310 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1957 103305283 : case t_FRAC:
1958 103305283 : if (!signe(x)) return gen_0;
1959 86530480 : z=cgetg(3,t_FRAC);
1960 86530667 : p1 = gcdii(x,gel(y,2));
1961 86529932 : if (equali1(p1))
1962 : {
1963 46541134 : set_avma((pari_sp)z);
1964 46541126 : gel(z,2) = icopy(gel(y,2));
1965 46541196 : gel(z,1) = mulii(gel(y,1), x);
1966 : }
1967 : else
1968 : {
1969 39989158 : x = diviiexact(x,p1); tetpil = avma;
1970 39988613 : gel(z,2) = diviiexact(gel(y,2), p1);
1971 39988361 : gel(z,1) = mulii(gel(y,1), x);
1972 39988833 : fix_frac_if_int_GC(z,tetpil);
1973 : }
1974 86530365 : return z;
1975 439334081 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1976 4199394 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1977 1904 : case t_QUAD: return mulRq(x,y);
1978 1607889 : case t_FFELT: return FF_Z_mul(y,x);
1979 : }
1980 :
1981 : case t_REAL:
1982 123857801 : switch(ty)
1983 : {
1984 38167518 : case t_FRAC: return mulrfrac(x, y);
1985 85690241 : case t_COMPLEX: return mulRc(x, y);
1986 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1987 21 : default: pari_err_TYPE2("*",x,y);
1988 : }
1989 :
1990 8288 : case t_INTMOD:
1991 8288 : switch(ty)
1992 : {
1993 7133 : case t_FRAC: { GEN X = gel(x,1);
1994 7133 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1995 7133 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1996 : }
1997 49 : case t_COMPLEX: return mulRc(x,y);
1998 511 : case t_PADIC: { GEN X = gel(x,1);
1999 511 : z = cgetg(3, t_INTMOD);
2000 511 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2001 : }
2002 63 : case t_QUAD: return mulRq(x, y);
2003 532 : case t_FFELT:
2004 532 : if (!equalii(gel(x,1),FF_p_i(y)))
2005 0 : pari_err_OP("*",x,y);
2006 532 : return FF_Z_mul(y,gel(x,2));
2007 : }
2008 :
2009 : case t_FRAC:
2010 7979431 : switch(ty)
2011 : {
2012 5394256 : case t_COMPLEX: return mulRc(x, y);
2013 2582752 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
2014 343 : case t_QUAD: return mulRq(x, y);
2015 2079 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
2016 : }
2017 :
2018 : case t_FFELT:
2019 36 : pari_err_TYPE2("*",x,y);
2020 :
2021 21 : case t_COMPLEX:
2022 21 : switch(ty)
2023 : {
2024 14 : case t_PADIC:
2025 14 : return Zp_nosquare_m1(padic_p(y))? mulRc(y, x): mulTp(x, y);
2026 7 : case t_QUAD:
2027 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
2028 7 : return mulqf(y, x, lx);
2029 : }
2030 :
2031 : case t_PADIC: /* ty == t_QUAD */
2032 28 : return (kro_quad(y,padic_p(x))== -1)? mulRq(x, y): mulTp(y, x);
2033 : }
2034 : }
2035 :
2036 166060461 : if (is_matvec_t(ty)) switch(tx)
2037 : {
2038 198686 : case t_VEC:
2039 198686 : switch(ty) {
2040 14854 : case t_COL: return RgV_RgC_mul(x,y);
2041 183832 : case t_MAT: return RgV_RgM_mul(x,y);
2042 : }
2043 0 : break;
2044 1687 : case t_COL:
2045 1687 : switch(ty) {
2046 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2047 0 : case t_MAT: return RgC_RgM_mul(x,y);
2048 : }
2049 0 : break;
2050 884061 : case t_MAT:
2051 884061 : switch(ty) {
2052 0 : case t_VEC: return RgM_RgV_mul(x,y);
2053 884061 : case t_COL: return RgM_RgC_mul(x,y);
2054 : }
2055 : default:
2056 5880383 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
2057 5880386 : return mul_self_scal(y, x);
2058 : }
2059 162616845 : if (is_matvec_t(tx))
2060 : {
2061 2146833 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2062 2147011 : return mul_self_scal(x, y);
2063 : }
2064 160470289 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2065 : /* tx < ty, !ismatvec(x and y) */
2066 :
2067 160470289 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2068 3265060 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2069 157205229 : if (is_scalar_t(tx)) {
2070 153737101 : if (tx == t_POLMOD) {
2071 3170919 : vx = varn(gel(x,1));
2072 3170919 : vy = gvar(y);
2073 3170919 : if (vx != vy) {
2074 2918576 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2075 23107 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2076 : }
2077 : /* error if ty == t_SER */
2078 252343 : av = avma; y = gmod(y, gel(x,1));
2079 252336 : return gc_upto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2080 : }
2081 150566182 : return mul_scal(y, x, ty);
2082 : }
2083 :
2084 : /* x and y are not scalars, nor matvec */
2085 3468114 : vx = gvar(x);
2086 3468126 : vy = gvar(y);
2087 3468126 : if (vx != vy) /* x or y is treated as a scalar */
2088 2785619 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2089 2785619 : : mul_scal(y, x, ty);
2090 : /* vx = vy */
2091 2044702 : switch(tx)
2092 : {
2093 2044667 : case t_POL:
2094 2044667 : switch (ty)
2095 : {
2096 6762 : case t_SER:
2097 : {
2098 : long e, E;
2099 6762 : av = avma; e = RgX_valrem(x, &x);
2100 6762 : if (e == LONG_MAX)
2101 21 : return gc_upto(av, gmul(Rg_get_0(x), Rg_get_0(y)));
2102 6741 : ly = lg(y); E = e + valser(y);
2103 6741 : if (ly == 2) { set_avma(av); return zeroser(vx, E); }
2104 6510 : if (ly == 3)
2105 : {
2106 : long s;
2107 : GEN c;
2108 63 : z = cgetg(3, t_SER); c = gmul(gel(x,2), gel(y,2));
2109 63 : s = !gequal0(c);
2110 63 : if (!s && signe(x) && signe(y) && isexactzero(c)) E++;
2111 63 : z[1] = evalsigne(s) | evalvarn(vx) | evalvalser(E);
2112 63 : gel(z,2) = c; return gc_upto(av, z);
2113 : }
2114 6447 : if (degpol(x))
2115 : {
2116 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2117 2030 : return gc_GEN(av, pol2ser(x, ly, E));
2118 : }
2119 : /* take advantage of x = c*t^e */
2120 4417 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2121 4417 : setvalser(y, e + valser(y)); return y;
2122 : }
2123 :
2124 2037905 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2125 : }
2126 0 : break;
2127 :
2128 35 : case t_SER:
2129 35 : switch (ty)
2130 : {
2131 35 : case t_RFRAC:
2132 35 : av = avma;
2133 35 : return gc_upto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2134 : }
2135 0 : break;
2136 : }
2137 0 : pari_err_TYPE2("*",x,y);
2138 : return NULL; /* LCOV_EXCL_LINE */
2139 : }
2140 :
2141 : /* return a nonnormalized result */
2142 : GEN
2143 114411 : sqr_ser_part(GEN x, long l1, long l2)
2144 : {
2145 : long i, j, l;
2146 : pari_sp av;
2147 : GEN Z, z, p1, p2;
2148 : long mi;
2149 114411 : if (l2 < l1) return zeroser(varn(x), 2*valser(x));
2150 114397 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2151 114397 : Z = cgetg(l2-l1+3, t_SER);
2152 114397 : Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
2153 114397 : z = Z + 2-l1;
2154 114397 : x += 2; mi = 0;
2155 427070 : for (i=0; i<l1; i++)
2156 : {
2157 312673 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2158 : }
2159 :
2160 725742 : for (i=l1; i<=l2; i++)
2161 : {
2162 611345 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2163 611345 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2164 1261895 : for (j=i-mi; j<=minss(l,mi); j++)
2165 650550 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2166 611345 : p1 = gshift(p1,1);
2167 611345 : if ((i&1) == 0 && p2[i>>1])
2168 95720 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2169 611345 : gel(z,i) = gc_upto(av,p1);
2170 : }
2171 114397 : return Z;
2172 : }
2173 :
2174 : /* (u + v X)^2 mod (X^2 + bX + c), b = 0 or -1
2175 : * = u^2 - c v^2 + (2uv - b v^2) X */
2176 : static GEN
2177 70 : sqrq(GEN x)
2178 : {
2179 70 : GEN T = gel(x,1), c = gel(T,2), b = gel(T,3);
2180 70 : GEN u = gel(x,2), v = gel(x,3), U, V, E, F, s, z;
2181 : pari_sp av, av2;
2182 :
2183 70 : z = cgetg(4, t_QUAD), gel(z,1) = ZX_copy(T); av = avma;
2184 70 : U = gsqr(u);
2185 70 : V = gsqr(v); s = gmul(c, V);
2186 70 : E = gmul(u, v);
2187 70 : F = signe(b)? gadd(E, V): E; /* u v - b v^2 */
2188 70 : av2 = avma;
2189 70 : gel(z,2) = gsub(U, s);
2190 70 : gel(z,3) = gadd(E, F);
2191 70 : gc_slice_unsafe(av, av2, z+2, 2); return z;
2192 : }
2193 :
2194 : GEN
2195 1217787972 : gsqr(GEN x)
2196 : {
2197 : long i, lx;
2198 : pari_sp av, tetpil;
2199 : GEN z, p1, p2, p3;
2200 :
2201 1217787972 : switch(typ(x))
2202 : {
2203 994721732 : case t_INT: return sqri(x);
2204 205545289 : case t_REAL: return sqrr(x);
2205 142425 : case t_INTMOD: { GEN X = gel(x,1);
2206 142425 : z = cgetg(3,t_INTMOD);
2207 142425 : gel(z,2) = gc_INT((pari_sp)z, remii(sqri(gel(x,2)), X));
2208 142426 : gel(z,1) = icopy(X); return z;
2209 : }
2210 4398046 : case t_FRAC: return sqrfrac(x);
2211 :
2212 8106610 : case t_COMPLEX:
2213 8106610 : if (isintzero(gel(x,1))) {
2214 586907 : av = avma;
2215 586907 : return gc_upto(av, gneg(gsqr(gel(x,2))));
2216 : }
2217 7519702 : z = cgetg(3,t_COMPLEX); av = avma;
2218 7519667 : p1 = gadd(gel(x,1),gel(x,2));
2219 7519428 : p2 = gsub(gel(x,1), gel(x,2));
2220 7519418 : p3 = gmul(gel(x,1),gel(x,2));
2221 7519575 : tetpil = avma;
2222 7519575 : gel(z,1) = gmul(p1,p2);
2223 7519593 : gel(z,2) = gshift(p3,1); gc_slice_unsafe(av,tetpil,z+1,2); return z;
2224 :
2225 4774 : case t_PADIC:
2226 : {
2227 4774 : GEN u = padic_u(x), p = padic_p(x), pd = padic_pd(x);
2228 4774 : long v2 = 2*valp(x), d = precp(x);
2229 4774 : if (!signe(u)) { x = gcopy(x); setvalp(x, v2); return x; }
2230 4718 : if (!absequaliu(p, 2))
2231 3185 : retmkpadic(Fp_sqr(u, pd), icopy(p), icopy(pd), v2, d);
2232 : /* p = 2*/
2233 1533 : if (d == 1) /* (1 + O(2))^2 = 1 + O(2^3) */
2234 7 : retmkpadic(gen_1, gen_2, utoipos(8), v2, 3);
2235 1526 : retmkpadic(remi2n(sqri(u), d + 1), gen_2, int2n(d + 1), v2, d + 1);
2236 : }
2237 70 : case t_QUAD: return sqrq(x);
2238 :
2239 426147 : case t_POLMOD: return sqr_polmod(gel(x,1), gel(x,2));
2240 :
2241 2971446 : case t_FFELT: return FF_sqr(x);
2242 :
2243 1461397 : case t_POL: return RgX_sqr(x);
2244 :
2245 36694 : case t_SER:
2246 36694 : lx = lg(x);
2247 36694 : if (ser_isexactzero(x)) {
2248 21 : GEN z = gcopy(x);
2249 21 : setvalser(z, 2*valser(x));
2250 21 : return z;
2251 : }
2252 36673 : if (lx < 40)
2253 36414 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2254 : else
2255 : {
2256 259 : pari_sp av = avma;
2257 259 : long v = 2 * valser(x);
2258 259 : x = ser2pol_i(x, lx);
2259 259 : x = RgXn_sqr(x, lx-2);
2260 259 : return gc_GEN(av, pol2ser(x, lx, v));
2261 : }
2262 :
2263 63 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2264 63 : gel(z,1) = gsqr(gel(x,1));
2265 63 : gel(z,2) = gsqr(gel(x,2)); return z;
2266 :
2267 1169 : case t_MAT: return RgM_sqr(x);
2268 28 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
2269 : /* fall through handle extended t_QFB */
2270 36 : case t_QFB: return qfbsqr(x);
2271 658 : case t_VECSMALL:
2272 658 : z = cgetg_copy(x, &lx);
2273 16289 : for (i=1; i<lx; i++)
2274 : {
2275 15631 : long xi = x[i];
2276 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2277 15631 : z[i] = x[xi];
2278 : }
2279 658 : return z;
2280 : }
2281 7 : pari_err_TYPE2("*",x,x);
2282 : return NULL; /* LCOV_EXCL_LINE */
2283 : }
2284 :
2285 : /********************************************************************/
2286 : /** **/
2287 : /** DIVISION **/
2288 : /** **/
2289 : /********************************************************************/
2290 : static GEN
2291 32139 : div_rfrac_scal(GEN x, GEN y)
2292 : {
2293 32139 : pari_sp av = avma;
2294 32139 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2295 32139 : return gc_upto(av, gred_rfrac_simple(gel(x,1), d));
2296 : }
2297 : static GEN
2298 37767 : div_scal_rfrac(GEN x, GEN y)
2299 : {
2300 37767 : GEN y1 = gel(y,1), y2 = gel(y,2);
2301 37767 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2302 : {
2303 14 : if (degpol(y1))
2304 : {
2305 14 : pari_sp av = avma;
2306 14 : GEN _1 = Rg_get_1(x);
2307 14 : if (base_change(_1)) y1 = gmul(y1, _1);
2308 14 : return gc_upto(av, gred_rfrac_simple(gmul(x, y2), y1));
2309 : }
2310 0 : y1 = gel(y1,2);
2311 : }
2312 37753 : return RgX_Rg_mul(y2, gdiv(x,y1));
2313 : }
2314 : static GEN
2315 1187251 : div_rfrac(GEN x, GEN y)
2316 1187251 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2317 :
2318 : /* x != 0 */
2319 : static GEN
2320 1608941 : div_ser_scal(GEN x, GEN t)
2321 : {
2322 1608941 : if (ser_isexactzero(x))
2323 : {
2324 28 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gdiv(gel(x,2), t), varn(x), 1);
2325 28 : setvalser(z, valser(x)); return z;
2326 : }
2327 7494879 : pari_APPLY_ser(gdiv(gel(x,i), t));
2328 : }
2329 : GEN
2330 7658 : ser_normalize(GEN x)
2331 : {
2332 7658 : long i, lx = lg(x);
2333 : GEN c, z;
2334 7658 : if (lx == 2) return x;
2335 7658 : c = gel(x,2); if (gequal1(c)) return x;
2336 7581 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2337 108752 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2338 7581 : return z;
2339 : }
2340 :
2341 : /* y != 0 */
2342 : static GEN
2343 7412314 : div_T_scal(GEN x, GEN y, long tx) {
2344 7412314 : switch(tx)
2345 : {
2346 5775084 : case t_POL: return RgX_Rg_div(x, y);
2347 1608934 : case t_SER: return div_ser_scal(x, y);
2348 28296 : case t_RFRAC: return div_rfrac_scal(x,y);
2349 : }
2350 0 : pari_err_TYPE2("/",x,y);
2351 : return NULL; /* LCOV_EXCL_LINE */
2352 : }
2353 :
2354 : static GEN
2355 9254006 : div_scal_pol(GEN x, GEN y) {
2356 9254006 : long ly = lg(y);
2357 : pari_sp av;
2358 : GEN _1;
2359 9254006 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2360 9180799 : if (isrationalzero(x)) return zeropol(varn(y));
2361 7133665 : av = avma;
2362 7133665 : _1 = Rg_get_1(x); if (base_change(_1)) y = gmul(y, _1);
2363 7133665 : return gc_upto(av, gred_rfrac_simple(x,y));
2364 : }
2365 : static GEN
2366 19292 : div_scal_ser(GEN x, GEN y)
2367 : {
2368 19292 : pari_sp av = avma;
2369 19292 : GEN _1 = Rg_get_1(x);
2370 19292 : if (base_change(_1)) y = gmul(y, _1);
2371 19292 : return gc_upto(av, gmul(x, ser_inv(y)));
2372 : }
2373 : static GEN
2374 9263698 : div_scal_T(GEN x, GEN y, long ty) {
2375 9263698 : switch(ty)
2376 : {
2377 9207402 : case t_POL: return div_scal_pol(x, y);
2378 19285 : case t_SER: return div_scal_ser(x, y);
2379 37011 : case t_RFRAC: return div_scal_rfrac(x, y);
2380 : }
2381 0 : pari_err_TYPE2("/",x,y);
2382 : return NULL; /* LCOV_EXCL_LINE */
2383 : }
2384 :
2385 : /* assume tx = ty = t_SER, same variable vx */
2386 : static GEN
2387 908987 : div_ser(GEN x, GEN y, long vx)
2388 : {
2389 908987 : long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
2390 908987 : GEN y0 = y, z;
2391 908987 : pari_sp av = avma;
2392 :
2393 908987 : if (!signe(y)) pari_err_INV("div_ser", y);
2394 908980 : if (ser_isexactzero(x))
2395 : {
2396 59892 : if (lx == 2) return zeroser(vx, v);
2397 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2398 147 : setvalser(z, v); return z;
2399 : }
2400 849088 : if (lx < ly) ly = lx;
2401 849088 : y = ser2pol_i_normalize(y, ly, &e);
2402 849088 : if (e)
2403 : {
2404 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2405 7 : ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
2406 : }
2407 849088 : if (ly == 3)
2408 : {
2409 12452 : z = cgetg(3, t_SER);
2410 12452 : z[1] = evalsigne(1) | evalvalser(v) | evalvarn(vx);
2411 12452 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2412 12452 : if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
2413 12452 : return gc_upto(av, z);
2414 : }
2415 836636 : x = ser2pol_i(x, ly);
2416 836636 : y = RgXn_div_i(x, y, ly-2);
2417 836636 : return gc_GEN(av, pol2ser(y, ly, v));
2418 : }
2419 : /* x,y compatible PADIC */
2420 : static GEN
2421 13212419 : divpp(GEN x, GEN y)
2422 : {
2423 13212419 : GEN M, ux = padic_u(x), uy = padic_u(y), p = padic_p(x);
2424 13212419 : long a = precp(x), b = precp(y), v = valp(x) - valp(y);
2425 :
2426 13212419 : if (!signe(uy)) pari_err_INV("divpp",y);
2427 13212419 : if (!signe(ux)) return zeropadic(p, v);
2428 13212090 : if (a > b) M = padic_pd(y); else { M = padic_pd(x); b = a; }
2429 13212090 : retmkpadic(Fp_div(ux, uy, M), icopy(padic_p(x)), icopy(M), v, b);
2430 : }
2431 : static GEN
2432 38087 : div_polmod_same(GEN T, GEN x, GEN y)
2433 : {
2434 38087 : long v = varn(T);
2435 38087 : GEN a, z = cgetg(3, t_POLMOD);
2436 38087 : gel(z,1) = RgX_copy(T);
2437 38087 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2438 12341 : a = gdiv(x, y);
2439 25746 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2440 119 : {
2441 119 : pari_sp av = avma;
2442 119 : a = gc_upto(av, gmul(x, RgXQ_inv(y, T)));
2443 : }
2444 25627 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2445 5348 : {
2446 5348 : pari_sp av = avma;
2447 5348 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2448 5348 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2449 5348 : a = gc_upto(av, a);
2450 : }
2451 : else
2452 : {
2453 20279 : pari_sp av = avma;
2454 20279 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2455 20279 : a = gc_upto(av, a);
2456 : }
2457 38087 : gel(z,2) = a; return z;
2458 : }
2459 : GEN
2460 464150565 : gdiv(GEN x, GEN y)
2461 : {
2462 464150565 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2463 : pari_sp av, tetpil;
2464 : GEN z, p1;
2465 :
2466 464150565 : if (tx == ty) switch(tx)
2467 : {
2468 94285387 : case t_INT:
2469 94285387 : return Qdivii(x,y);
2470 :
2471 109130224 : case t_REAL: return divrr(x,y);
2472 25261 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2473 25261 : z = cgetg(3,t_INTMOD);
2474 25261 : if (X==Y || equalii(X,Y))
2475 25247 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2476 14 : gel(z,1) = gcdii(X,Y);
2477 14 : warn_coercion(X,Y,gel(z,1));
2478 14 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2479 14 : gel(z,2) = gc_INT(av, remii(p1, gel(z,1))); return z;
2480 : }
2481 4212459 : case t_FRAC: {
2482 4212459 : GEN x1 = gel(x,1), x2 = gel(x,2);
2483 4212459 : GEN y1 = gel(y,1), y2 = gel(y,2);
2484 4212459 : z = cgetg(3, t_FRAC);
2485 4212458 : p1 = gcdii(x1, y1);
2486 4212459 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2487 4212459 : p1 = gcdii(x2, y2);
2488 4212457 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2489 4212450 : tetpil = avma;
2490 4212450 : gel(z,2) = mulii(x2,y1);
2491 4212452 : gel(z,1) = mulii(x1,y2);
2492 4212454 : normalize_frac(z);
2493 4212453 : fix_frac_if_int_GC(z,tetpil);
2494 4212458 : return z;
2495 : }
2496 10696426 : case t_COMPLEX:
2497 10696426 : if (isintzero(gel(y,1)))
2498 : {
2499 486220 : y = gel(y,2);
2500 486220 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2501 277424 : z = cgetg(3,t_COMPLEX);
2502 277424 : gel(z,1) = gdiv(gel(x,2), y);
2503 277424 : av = avma;
2504 277424 : gel(z,2) = gc_upto(av, gneg(gdiv(gel(x,1), y)));
2505 277424 : return z;
2506 : }
2507 10210204 : av = avma;
2508 10210204 : return gc_upto(av, gdiv(mulcc(x, conj_i(y)), cxnorm(y)));
2509 :
2510 587462 : case t_PADIC:
2511 587462 : if (!equalii(padic_p(x), padic_p(y))) pari_err_OP("/",x,y);
2512 587462 : return divpp(x, y);
2513 :
2514 259 : case t_QUAD:
2515 259 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2516 245 : av = avma;
2517 245 : return gc_upto(av, gdiv(mulqq(x, conj_i(y)), quadnorm(y)));
2518 :
2519 2464024 : case t_FFELT: return FF_div(x,y);
2520 :
2521 42462 : case t_POLMOD:
2522 42462 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2523 38087 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2524 : else {
2525 4375 : av = avma;
2526 4375 : z = gc_upto(av, gmul(x, ginv(y)));
2527 : }
2528 42462 : return z;
2529 :
2530 18567558 : case t_POL:
2531 18567558 : vx = varn(x);
2532 18567558 : vy = varn(y);
2533 18567558 : if (vx != vy) {
2534 154084 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2535 46604 : else return div_scal_pol(x, y);
2536 : }
2537 18413474 : if (!signe(y)) pari_err_INV("gdiv",y);
2538 18413474 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2539 18288725 : av = avma;
2540 18288725 : return gc_upto(av, gred_rfrac(x,y));
2541 :
2542 909008 : case t_SER:
2543 909008 : vx = varn(x);
2544 909008 : vy = varn(y);
2545 909008 : if (vx != vy) {
2546 21 : if (varncmp(vx, vy) < 0)
2547 : {
2548 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2549 7 : return div_ser_scal(x, y);
2550 : }
2551 7 : return div_scal_ser(x, y);
2552 : }
2553 908987 : return div_ser(x, y, vx);
2554 1191710 : case t_RFRAC:
2555 1191710 : vx = varn(gel(x,2));
2556 1191710 : vy = varn(gel(y,2));
2557 1191710 : if (vx != vy) {
2558 4459 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2559 756 : else return div_scal_rfrac(x, y);
2560 : }
2561 1187251 : return div_rfrac(x,y);
2562 :
2563 21 : case t_VEC: /* handle extended t_QFB */
2564 21 : case t_QFB: av = avma; return gc_upto(av, qfbcomp(x, ginv(y)));
2565 :
2566 28 : case t_MAT:
2567 28 : p1 = RgM_div(x,y);
2568 28 : if (!p1) pari_err_INV("gdiv",y);
2569 21 : return p1;
2570 :
2571 0 : default: pari_err_TYPE2("/",x,y);
2572 : }
2573 :
2574 222042300 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2575 : {
2576 3937325 : long s = signe(x);
2577 3937325 : if (!s) {
2578 529904 : if (gequal0(y)) pari_err_INV("gdiv",y);
2579 529904 : switch (ty)
2580 : {
2581 526761 : default: return gen_0;
2582 28 : case t_INTMOD:
2583 28 : z = cgetg(3,t_INTMOD);
2584 28 : gel(z,1) = icopy(gel(y,1));
2585 28 : gel(z,2) = gen_0; return z;
2586 3115 : case t_FFELT: return FF_zero(y);
2587 : }
2588 : }
2589 3407421 : if (is_pm1(x)) {
2590 1328283 : if (s > 0) return ginv(y);
2591 230708 : av = avma; return gc_upto(av, ginv(gneg(y)));
2592 : }
2593 2079137 : switch(ty)
2594 : {
2595 663756 : case t_REAL: return divir(x,y);
2596 21 : case t_INTMOD:
2597 21 : z = cgetg(3, t_INTMOD);
2598 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2599 875171 : case t_FRAC:
2600 875171 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2601 875171 : if (equali1(p1))
2602 : {
2603 354182 : set_avma((pari_sp)z);
2604 354182 : gel(z,2) = icopy(gel(y,1));
2605 354182 : gel(z,1) = mulii(gel(y,2), x);
2606 354182 : normalize_frac(z);
2607 354182 : fix_frac_if_int(z);
2608 : }
2609 : else
2610 : {
2611 520989 : x = diviiexact(x,p1); tetpil = avma;
2612 520989 : gel(z,2) = diviiexact(gel(y,1), p1);
2613 520989 : gel(z,1) = mulii(gel(y,2), x);
2614 520989 : normalize_frac(z);
2615 520989 : fix_frac_if_int_GC(z,tetpil);
2616 : }
2617 875171 : return z;
2618 :
2619 427 : case t_FFELT: return Z_FF_div(x,y);
2620 538032 : case t_COMPLEX: return divRc(x,y);
2621 1505 : case t_PADIC: return divTp(x, y);
2622 231 : case t_QUAD: return divRq(x,y);
2623 : }
2624 : }
2625 218104969 : if (gequal0(y))
2626 : {
2627 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2628 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2629 : }
2630 :
2631 218113685 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2632 : {
2633 175805650 : case t_REAL:
2634 175805650 : switch(ty)
2635 : {
2636 172417895 : case t_INT: return divri(x,y);
2637 518709 : case t_FRAC:
2638 518709 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2639 518709 : return gc_leaf(av, z);
2640 2869112 : case t_COMPLEX: return divRc(x, y);
2641 42 : case t_QUAD: return divfq(x, y, realprec(x));
2642 18 : default: pari_err_TYPE2("/",x,y);
2643 : }
2644 :
2645 301 : case t_INTMOD:
2646 301 : switch(ty)
2647 : {
2648 217 : case t_INT:
2649 217 : z = cgetg(3, t_INTMOD);
2650 217 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2651 28 : case t_FRAC: { GEN X = gel(x,1);
2652 28 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2653 28 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2654 : }
2655 7 : case t_FFELT:
2656 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2657 0 : pari_err_OP("/",x,y);
2658 7 : return Z_FF_div(gel(x,2),y);
2659 :
2660 0 : case t_COMPLEX: return divRc(x,y);
2661 14 : case t_QUAD: return divRq(x,y);
2662 :
2663 7 : case t_PADIC: { GEN X = gel(x,1);
2664 7 : z = cgetg(3, t_INTMOD);
2665 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2666 : }
2667 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2668 : }
2669 :
2670 : case t_FRAC:
2671 2186499 : switch(ty)
2672 : {
2673 2106355 : case t_INT: z = cgetg(3, t_FRAC);
2674 2106355 : p1 = gcdii(y,gel(x,1));
2675 2106355 : if (equali1(p1))
2676 : {
2677 837348 : set_avma((pari_sp)z); tetpil = 0;
2678 837348 : gel(z,1) = icopy(gel(x,1));
2679 : }
2680 : else
2681 : {
2682 1269007 : y = diviiexact(y,p1); tetpil = avma;
2683 1269007 : gel(z,1) = diviiexact(gel(x,1), p1);
2684 : }
2685 2106355 : gel(z,2) = mulii(gel(x,2),y);
2686 2106355 : normalize_frac(z);
2687 2106355 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2688 2106355 : return z;
2689 :
2690 77058 : case t_REAL:
2691 77058 : av = avma;
2692 77058 : return gc_leaf(av, divir(gel(x,1), mulri(y,gel(x,2))));
2693 :
2694 7 : case t_INTMOD: { GEN Y = gel(y,1);
2695 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2696 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2697 : }
2698 :
2699 28 : case t_FFELT: av=avma;
2700 28 : return gc_upto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2701 :
2702 868 : case t_COMPLEX: return divRc(x, y);
2703 :
2704 2141 : case t_PADIC:
2705 2141 : if (!signe(gel(x,1))) return gen_0;
2706 2141 : return divTp(x, y);
2707 :
2708 42 : case t_QUAD: return divRq(x, y);
2709 : }
2710 :
2711 : case t_FFELT:
2712 161 : switch (ty)
2713 : {
2714 77 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2715 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2716 7 : case t_INTMOD:
2717 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2718 0 : pari_err_OP("/",x,y);
2719 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2720 49 : default:
2721 49 : pari_err_TYPE2("/",x,y);
2722 : }
2723 0 : break;
2724 :
2725 15064603 : case t_COMPLEX:
2726 15064603 : switch(ty)
2727 : {
2728 15064605 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2729 0 : case t_INTMOD: return mulRc(ginv(y), x);
2730 0 : case t_PADIC:
2731 0 : return Zp_nosquare_m1(padic_p(y))? divcR(x,y): divTp(x, y);
2732 0 : case t_QUAD:
2733 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2734 0 : return divfq(x, y, lx);
2735 : }
2736 :
2737 : case t_PADIC:
2738 1205055 : switch(ty)
2739 : {
2740 1205001 : case t_INT: case t_FRAC: { GEN p = padic_p(x);
2741 1204910 : return signe(padic_u(x))? divpT(x, y)
2742 2409911 : : zeropadic(p, valp(x) - Q_pval(y,p));
2743 : }
2744 7 : case t_INTMOD: { GEN Y = gel(y,1);
2745 7 : z = cgetg(3, t_INTMOD);
2746 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2747 : }
2748 0 : case t_COMPLEX: return divRc(x,y);
2749 14 : case t_QUAD: return divRq(x,y);
2750 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2751 : }
2752 :
2753 : case t_QUAD:
2754 1545 : switch (ty)
2755 : {
2756 1435 : case t_INT: case t_INTMOD: case t_FRAC:
2757 1435 : z = cgetg(4,t_QUAD);
2758 1435 : gel(z,1) = ZX_copy(gel(x,1));
2759 1435 : gel(z,2) = gdiv(gel(x,2), y);
2760 1435 : gel(z,3) = gdiv(gel(x,3), y); return z;
2761 63 : case t_REAL: return divqf(x, y, realprec(y));
2762 28 : case t_PADIC: return divTp(x, y);
2763 0 : case t_COMPLEX:
2764 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2765 0 : return divqf(x, y, ly);
2766 : }
2767 : }
2768 23849789 : switch(ty) {
2769 588777 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2770 588777 : av = avma; return gc_upto(av, gmul(x, ginv(y)));
2771 42 : case t_MAT:
2772 42 : av = avma; p1 = RgM_inv(y);
2773 35 : if (!p1) pari_err_INV("gdiv",y);
2774 28 : return gc_upto(av, gmul(x, p1));
2775 0 : case t_VEC: case t_COL:
2776 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2777 0 : pari_err_TYPE2("/",x,y);
2778 : }
2779 23260970 : switch(tx) {
2780 4427591 : case t_VEC: case t_COL: case t_MAT:
2781 14382595 : pari_APPLY_same(gdiv(gel(x,i),y));
2782 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2783 0 : pari_err_TYPE2("/",x,y);
2784 : }
2785 :
2786 18833379 : vy = gvar(y);
2787 18834352 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2788 218434 : vx = varn(X);
2789 218434 : if (vx != vy) {
2790 217664 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2791 217251 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2792 : }
2793 : /* y is POL, SER or RFRAC */
2794 770 : av = avma;
2795 770 : switch(ty)
2796 : {
2797 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2798 770 : default: y = ginvmod(gmod(y,X), X);
2799 : }
2800 763 : return gc_upto(av, mul_polmod_same(X, gel(x,2), y));
2801 : }
2802 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2803 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2804 : * variable NO_VARIABLE, then the operation is incorrect */
2805 18615918 : vx = gvar(x);
2806 18615905 : if (vx != vy) { /* includes cases where one is scalar */
2807 16675607 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2808 9263285 : else return div_scal_T(x, y, ty);
2809 : }
2810 1940298 : switch(tx)
2811 : {
2812 1046185 : case t_POL:
2813 1046185 : switch(ty)
2814 : {
2815 28 : case t_SER:
2816 : {
2817 28 : GEN y0 = y;
2818 : long v;
2819 28 : av = avma; v = RgX_valrem(x, &x);
2820 28 : if (v == LONG_MAX) return gc_upto(av, Rg_get_0(x));
2821 14 : v -= valser(y); ly = lg(y); /* > 2 */
2822 14 : y = ser2pol_i_normalize(y, ly, &i);
2823 14 : if (i)
2824 : {
2825 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2826 7 : ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
2827 : }
2828 7 : x = RgXn_div(x, y, ly-2);
2829 7 : return gc_GEN(av, pol2ser(x, ly, v));
2830 : }
2831 :
2832 1046157 : case t_RFRAC:
2833 : {
2834 1046157 : GEN y1 = gel(y,1), y2 = gel(y,2);
2835 1046157 : if (typ(y1) == t_POL && varn(y1) == vx)
2836 1157 : return mul_rfrac_scal(y2, y1, x);
2837 1045000 : av = avma;
2838 1045000 : return gc_upto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2839 : }
2840 : }
2841 0 : break;
2842 :
2843 631728 : case t_SER:
2844 631728 : switch(ty)
2845 : {
2846 631714 : case t_POL:
2847 : {
2848 631714 : long e = valser(x);
2849 631714 : lx = lg(x);
2850 631714 : if (lx == 2) return zeroser(vx, e - RgX_val(y));
2851 631609 : av = avma; e -= RgX_valrem_inexact(y, &y);
2852 631609 : if (lx == 3)
2853 : {
2854 : GEN c;
2855 51943 : z = cgetg(3, t_SER); c = gdiv(gel(x,2), gel(y,2));
2856 51943 : z[1] = evalsigne(!gequal0(c)) | evalvarn(vx) | evalvalser(e);
2857 51943 : gel(z,2) = c; return z;
2858 : }
2859 579666 : x = ser2pol_i(x, lx);
2860 579666 : x = RgXn_div(x, y, lx - 2);
2861 579666 : return gc_GEN(av, pol2ser(x, lx, e));
2862 : }
2863 14 : case t_RFRAC:
2864 14 : av = avma;
2865 14 : return gc_upto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2866 : }
2867 0 : break;
2868 :
2869 262365 : case t_RFRAC:
2870 262365 : switch(ty)
2871 : {
2872 262365 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2873 0 : case t_SER:
2874 0 : av = avma;
2875 0 : return gc_upto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2876 : }
2877 0 : break;
2878 : }
2879 20 : pari_err_TYPE2("/",x,y);
2880 : return NULL; /* LCOV_EXCL_LINE */
2881 : }
2882 :
2883 : /********************************************************************/
2884 : /** **/
2885 : /** SIMPLE MULTIPLICATION **/
2886 : /** **/
2887 : /********************************************************************/
2888 : GEN
2889 261418049 : gmulsg(long s, GEN x)
2890 : {
2891 : pari_sp av;
2892 : long i;
2893 : GEN z;
2894 :
2895 261418049 : switch(typ(x))
2896 : {
2897 146912911 : case t_INT: return mulsi(s,x);
2898 83225146 : case t_REAL: return s? mulsr(s,x): gen_0; /* gmul semantic */
2899 174800 : case t_INTMOD: { GEN p = gel(x,1);
2900 174800 : z = cgetg(3,t_INTMOD);
2901 174800 : gel(z,2) = gc_INT((pari_sp)z, modii(mulsi(s,gel(x,2)), p));
2902 174799 : gel(z,1) = icopy(p); return z;
2903 : }
2904 2337937 : case t_FFELT: return FF_Z_mul(x,stoi(s));
2905 6163363 : case t_FRAC:
2906 6163363 : if (!s) return gen_0;
2907 6152639 : z = cgetg(3,t_FRAC);
2908 6152639 : i = labs(s); i = ugcdiu(gel(x,2), i);
2909 6152639 : if (i == 1)
2910 : {
2911 2263664 : gel(z,2) = icopy(gel(x,2));
2912 2263664 : gel(z,1) = mulis(gel(x,1), s);
2913 : }
2914 : else
2915 : {
2916 3888975 : gel(z,2) = diviuexact(gel(x,2), (ulong)i);
2917 3888975 : gel(z,1) = mulis(gel(x,1), s/i);
2918 3888975 : fix_frac_if_int(z);
2919 : }
2920 6152639 : return z;
2921 :
2922 20410183 : case t_COMPLEX:
2923 20410183 : if (!s) return gen_0;
2924 18551679 : z = cgetg(3, t_COMPLEX);
2925 18552487 : gel(z,1) = gmulsg(s,gel(x,1));
2926 18580646 : gel(z,2) = gmulsg(s,gel(x,2)); return z;
2927 :
2928 1449 : case t_PADIC:
2929 1449 : if (!s) return gen_0;
2930 1449 : av = avma; return gc_upto(av, mulpp(cvtop2(stoi(s),x), x));
2931 :
2932 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2933 7 : gel(z,1) = ZX_copy(gel(x,1));
2934 7 : gel(z,2) = gmulsg(s,gel(x,2));
2935 7 : gel(z,3) = gmulsg(s,gel(x,3)); return z;
2936 :
2937 220724 : case t_POLMOD:
2938 220724 : retmkpolmod(gmulsg(s,gel(x,2)), RgX_copy(gel(x,1)));
2939 :
2940 827309 : case t_POL:
2941 827309 : if (!signe(x)) return RgX_copy(x);
2942 808332 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2943 3339699 : pari_APPLY_pol(gmulsg(s,gel(x,i)));
2944 :
2945 182 : case t_SER:
2946 182 : if (ser_isexactzero(x)) return gcopy(x);
2947 182 : if (!s) return Rg_get_0(x);
2948 3864 : pari_APPLY_ser(gmulsg(s,gel(x,i)));
2949 :
2950 0 : case t_RFRAC:
2951 0 : if (!s) return zeropol(varn(gel(x,2)));
2952 0 : if (s == 1) return gcopy(x);
2953 0 : if (s == -1) return gneg(x);
2954 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), stoi(s));
2955 :
2956 1306368 : case t_VEC: case t_COL: case t_MAT:
2957 4116772 : pari_APPLY_same(gmulsg(s,gel(x,i)));
2958 : }
2959 0 : pari_err_TYPE("gmulsg",x);
2960 : return NULL; /* LCOV_EXCL_LINE */
2961 : }
2962 :
2963 : GEN
2964 283956448 : gmulug(ulong s, GEN x)
2965 : {
2966 : pari_sp av;
2967 : long i;
2968 : GEN z;
2969 :
2970 283956448 : switch(typ(x))
2971 : {
2972 281555476 : case t_INT: return mului(s,x);
2973 1068534 : case t_REAL: return s? mulur(s,x): gen_0; /* gmul semantic */
2974 364 : case t_INTMOD: { GEN p = gel(x,1);
2975 364 : z = cgetg(3,t_INTMOD);
2976 364 : gel(z,2) = gc_INT((pari_sp)z, modii(mului(s,gel(x,2)), p));
2977 364 : gel(z,1) = icopy(p); return z;
2978 : }
2979 413 : case t_FFELT: return FF_Z_mul(x,utoi(s));
2980 1064129 : case t_FRAC:
2981 1064129 : if (!s) return gen_0;
2982 1064115 : z = cgetg(3,t_FRAC);
2983 1064115 : i = ugcdiu(gel(x,2), s);
2984 1064115 : if (i == 1)
2985 : {
2986 905103 : gel(z,2) = icopy(gel(x,2));
2987 905103 : gel(z,1) = muliu(gel(x,1), s);
2988 : }
2989 : else
2990 : {
2991 159012 : gel(z,2) = diviuexact(gel(x,2), i);
2992 159012 : gel(z,1) = muliu(gel(x,1), s/i);
2993 159012 : fix_frac_if_int(z);
2994 : }
2995 1064115 : return z;
2996 :
2997 35973 : case t_COMPLEX:
2998 35973 : if (!s) return gen_0;
2999 35973 : z = cgetg(3, t_COMPLEX);
3000 35973 : gel(z,1) = gmulug(s,gel(x,1));
3001 35973 : gel(z,2) = gmulug(s,gel(x,2)); return z;
3002 :
3003 7342 : case t_PADIC:
3004 7342 : if (!s) return gen_0;
3005 7342 : av = avma; return gc_upto(av, mulpp(cvtop2(utoi(s),x), x));
3006 :
3007 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3008 0 : gel(z,1) = ZX_copy(gel(x,1));
3009 0 : gel(z,2) = gmulug(s,gel(x,2));
3010 0 : gel(z,3) = gmulug(s,gel(x,3)); return z;
3011 :
3012 6363 : case t_POLMOD:
3013 6363 : retmkpolmod(gmulug(s,gel(x,2)), RgX_copy(gel(x,1)));
3014 :
3015 200235 : case t_POL:
3016 200235 : if (!signe(x)) return RgX_copy(x);
3017 199920 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
3018 621229 : pari_APPLY_pol(gmulug(s,gel(x,i)));
3019 :
3020 0 : case t_SER:
3021 0 : if (ser_isexactzero(x)) return gcopy(x);
3022 0 : if (!s) return Rg_get_0(x);
3023 0 : pari_APPLY_ser(gmulug(s,gel(x,i)));
3024 :
3025 0 : case t_RFRAC:
3026 0 : if (!s) return zeropol(varn(gel(x,2)));
3027 0 : if (s == 1) return gcopy(x);
3028 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), utoi(s));
3029 :
3030 17619 : case t_VEC: case t_COL: case t_MAT:
3031 96495 : pari_APPLY_same(gmulug(s,gel(x,i)));
3032 : }
3033 0 : pari_err_TYPE("gmulsg",x);
3034 : return NULL; /* LCOV_EXCL_LINE */
3035 : }
3036 :
3037 : /********************************************************************/
3038 : /** **/
3039 : /** SIMPLE DIVISION **/
3040 : /** **/
3041 : /********************************************************************/
3042 :
3043 : GEN
3044 29176177 : gdivgs(GEN x, long s)
3045 : {
3046 29176177 : long tx = typ(x), i;
3047 : pari_sp av;
3048 : GEN z;
3049 :
3050 29176177 : if (!s)
3051 : {
3052 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3053 0 : pari_err_INV("gdivgs",gen_0);
3054 : }
3055 29321813 : switch(tx)
3056 : {
3057 2582173 : case t_INT: return Qdivis(x, s);
3058 18729015 : case t_REAL: return divrs(x,s);
3059 :
3060 357 : case t_INTMOD:
3061 357 : z = cgetg(3, t_INTMOD);
3062 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3063 :
3064 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3065 :
3066 641541 : case t_FRAC: z = cgetg(3, t_FRAC);
3067 641541 : i = ugcdiu(gel(x,1), labs(s));
3068 641541 : if (i == 1)
3069 : {
3070 482891 : gel(z,2) = mulsi(s, gel(x,2));
3071 482891 : gel(z,1) = icopy(gel(x,1));
3072 : }
3073 : else
3074 : {
3075 158650 : gel(z,2) = mulsi(s/i, gel(x,2));
3076 158650 : gel(z,1) = divis(gel(x,1), i);
3077 : }
3078 641541 : normalize_frac(z);
3079 641541 : fix_frac_if_int(z); return z;
3080 :
3081 6884903 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3082 6866596 : gel(z,1) = gdivgs(gel(x,1),s);
3083 6893773 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3084 :
3085 133 : case t_PADIC: /* divpT */
3086 : {
3087 133 : GEN p = padic_p(x);
3088 133 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3089 133 : av = avma;
3090 133 : return gc_upto(av, divpp(x, cvtop2(stoi(s),x)));
3091 : }
3092 :
3093 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3094 28 : gel(z,1) = ZX_copy(gel(x,1));
3095 28 : gel(z,2) = gdivgs(gel(x,2),s);
3096 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3097 :
3098 37282 : case t_POLMOD:
3099 37282 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3100 :
3101 91 : case t_RFRAC:
3102 91 : if (s == 1) return gcopy(x);
3103 84 : else if (s == -1) return gneg(x);
3104 84 : return div_rfrac_scal(x, stoi(s));
3105 :
3106 156492 : case t_POL: pari_APPLY_pol_normalized(gdivgs(gel(x,i),s));
3107 0 : case t_SER: pari_APPLY_ser_normalized(gdivgs(gel(x,i),s));
3108 408630 : case t_VEC:
3109 : case t_COL:
3110 897428 : case t_MAT: pari_APPLY_same(gdivgs(gel(x,i),s));
3111 : }
3112 0 : pari_err_TYPE2("/",x, stoi(s));
3113 : return NULL; /* LCOV_EXCL_LINE */
3114 : }
3115 :
3116 : GEN
3117 61977051 : gdivgu(GEN x, ulong s)
3118 : {
3119 61977051 : long tx = typ(x), i;
3120 : pari_sp av;
3121 : GEN z;
3122 :
3123 61977051 : if (!s)
3124 : {
3125 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3126 0 : pari_err_INV("gdivgu",gen_0);
3127 : }
3128 61976924 : switch(tx)
3129 : {
3130 25111598 : case t_INT: return Qdiviu(x, s);
3131 12923036 : case t_REAL: return divru(x,s);
3132 :
3133 210315 : case t_INTMOD:
3134 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3135 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3136 :
3137 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3138 :
3139 728375 : case t_FRAC: z = cgetg(3, t_FRAC);
3140 728375 : i = ugcdiu(gel(x,1), s);
3141 728375 : if (i == 1)
3142 : {
3143 535476 : gel(z,2) = mului(s, gel(x,2));
3144 535476 : gel(z,1) = icopy(gel(x,1));
3145 : }
3146 : else
3147 : {
3148 192899 : gel(z,2) = mului(s/i, gel(x,2));
3149 192899 : gel(z,1) = divis(gel(x,1), i);
3150 : }
3151 728375 : normalize_frac(z);
3152 728375 : fix_frac_if_int(z); return z;
3153 :
3154 11347580 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3155 11347580 : gel(z,1) = gdivgu(gel(x,1),s);
3156 11347581 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3157 :
3158 11551006 : case t_PADIC: /* divpT */
3159 : {
3160 11551006 : GEN p = padic_p(x);
3161 11551006 : if (!signe(padic_u(x))) return zeropadic(p, valp(x) - u_pval(s,p));
3162 11415447 : av = avma;
3163 11415447 : return gc_upto(av, divpp(x, cvtop2(utoi(s),x)));
3164 : }
3165 :
3166 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3167 0 : gel(z,1) = ZX_copy(gel(x,1));
3168 0 : gel(z,2) = gdivgu(gel(x,2),s);
3169 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3170 :
3171 1456 : case t_POLMOD:
3172 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3173 :
3174 56 : case t_RFRAC:
3175 56 : if (s == 1) return gcopy(x);
3176 56 : return div_rfrac_scal(x, utoi(s));
3177 :
3178 472878 : case t_POL: pari_APPLY_pol_normalized(gdivgu(gel(x,i),s));
3179 18389 : case t_SER: pari_APPLY_ser_normalized(gdivgu(gel(x,i),s));
3180 336 : case t_VEC:
3181 : case t_COL:
3182 1162 : case t_MAT: pari_APPLY_same(gdivgu(gel(x,i),s));
3183 : }
3184 0 : pari_err_TYPE2("/",x, utoi(s));
3185 : return NULL; /* LCOV_EXCL_LINE */
3186 : }
3187 :
3188 : /* x / (i*(i+1)) */
3189 : GEN
3190 214278705 : divrunextu(GEN x, ulong i)
3191 : {
3192 214278705 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3193 0 : return divri(x, muluu(i, i+1));
3194 : else
3195 214278705 : return divru(x, i*(i+1));
3196 : }
3197 : /* x / (i*(i+1)) */
3198 : GEN
3199 1583943 : gdivgunextu(GEN x, ulong i)
3200 : {
3201 1583943 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3202 0 : return gdivgu(x, i*(i+1));
3203 : else
3204 1583943 : return gdiv(x, muluu(i, i+1));
3205 : }
3206 :
3207 : /* True shift (exact multiplication by 2^n) */
3208 : GEN
3209 155001957 : gmul2n(GEN x, long n)
3210 : {
3211 : GEN z, a, b;
3212 : long k, l;
3213 :
3214 155001957 : switch(typ(x))
3215 : {
3216 44010188 : case t_INT:
3217 44010188 : if (n>=0) return shifti(x,n);
3218 7896824 : if (!signe(x)) return gen_0;
3219 5206220 : l = vali(x); n = -n;
3220 5206281 : if (n<=l) return shifti(x,-n);
3221 396713 : z = cgetg(3,t_FRAC);
3222 396713 : gel(z,1) = shifti(x,-l);
3223 396713 : gel(z,2) = int2n(n-l); return z;
3224 :
3225 70278411 : case t_REAL:
3226 70278411 : return shiftr(x,n);
3227 :
3228 294809 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3229 294809 : z = cgetg(3,t_INTMOD);
3230 294809 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3231 196648 : gel(z,2) = gc_INT((pari_sp)z, modii(shifti(a,n), b));
3232 196645 : gel(z,1) = icopy(b); return z;
3233 :
3234 3843732 : case t_FFELT: return FF_mul2n(x,n);
3235 :
3236 1352609 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3237 1352609 : l = vali(a);
3238 1352609 : k = vali(b);
3239 1352609 : if (n+l >= k)
3240 : {
3241 472609 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3242 376258 : l = n-k; k = -k;
3243 : }
3244 : else
3245 : {
3246 880000 : k = -(l+n); l = -l;
3247 : }
3248 1256258 : z = cgetg(3,t_FRAC);
3249 1256258 : gel(z,1) = shifti(a,l);
3250 1256257 : gel(z,2) = shifti(b,k); return z;
3251 :
3252 11434693 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3253 11434687 : gel(z,1) = gmul2n(gel(x,1),n);
3254 11434682 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3255 :
3256 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3257 105 : gel(z,1) = ZX_copy(gel(x,1));
3258 105 : gel(z,2) = gmul2n(gel(x,2),n);
3259 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3260 :
3261 178227 : case t_POLMOD:
3262 178227 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3263 :
3264 6871001 : case t_POL:
3265 24796888 : pari_APPLY_pol(gmul2n(gel(x,i),n));
3266 106459 : case t_SER:
3267 106459 : if (ser_isexactzero(x)) return gcopy(x);
3268 650559 : pari_APPLY_ser(gmul2n(gel(x,i),n));
3269 16625999 : case t_VEC: case t_COL: case t_MAT:
3270 61609219 : pari_APPLY_same(gmul2n(gel(x,i),n));
3271 :
3272 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3273 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3274 :
3275 7644 : case t_PADIC: /* int2n wrong if n < 0 */
3276 7644 : return gmul(gmul2n(gen_1,n),x);
3277 : }
3278 0 : pari_err_TYPE("gmul2n",x);
3279 : return NULL; /* LCOV_EXCL_LINE */
3280 : }
3281 :
3282 : /*******************************************************************/
3283 : /* */
3284 : /* INVERSE */
3285 : /* */
3286 : /*******************************************************************/
3287 : static GEN
3288 192251 : inv_polmod(GEN T, GEN x)
3289 : {
3290 192251 : GEN z = cgetg(3,t_POLMOD), a;
3291 192250 : gel(z,1) = RgX_copy(T);
3292 192250 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3293 67288 : a = ginv(x);
3294 : else
3295 : {
3296 124962 : if (lg(T) == 5) /* quadratic fields */
3297 13377 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3298 : else
3299 111585 : a = RgXQ_inv(x, T);
3300 : }
3301 192250 : gel(z,2) = a; return z;
3302 : }
3303 :
3304 : GEN
3305 38787012 : ginv(GEN x)
3306 : {
3307 : long s;
3308 : pari_sp av;
3309 : GEN z, y;
3310 :
3311 38787012 : switch(typ(x))
3312 : {
3313 9718765 : case t_INT:
3314 9718765 : if (is_pm1(x)) return icopy(x);
3315 7998441 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3316 7998427 : z = cgetg(3,t_FRAC);
3317 8005248 : gel(z,1) = s<0? gen_m1: gen_1;
3318 8005248 : gel(z,2) = absi(x); return z;
3319 :
3320 12656533 : case t_REAL: return invr(x);
3321 :
3322 6517 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3323 6517 : gel(z,1) = icopy(gel(x,1));
3324 6517 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3325 :
3326 554095 : case t_FRAC: {
3327 554095 : GEN a = gel(x,1), b = gel(x,2);
3328 554095 : s = signe(a);
3329 554095 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3330 263141 : z = cgetg(3,t_FRAC);
3331 263141 : gel(z,1) = icopy(b);
3332 263144 : gel(z,2) = icopy(a);
3333 263145 : normalize_frac(z); return z;
3334 : }
3335 9736770 : case t_COMPLEX:
3336 9736770 : av = avma;
3337 9736770 : return gc_upto(av, divcR(conj_i(x), cxnorm(x)));
3338 :
3339 497 : case t_QUAD:
3340 497 : av = avma;
3341 497 : return gc_upto(av, gdiv(conj_i(x), quadnorm(x)));
3342 :
3343 358358 : case t_PADIC:
3344 : {
3345 358358 : GEN p = padic_p(x), pd = padic_pd(x), u = padic_u(x);
3346 358358 : long d = precp(x);
3347 358358 : if (!signe(u)) pari_err_INV("ginv",x);
3348 358351 : retmkpadic(Zp_inv(u, p, d), icopy(p), icopy(pd), -valp(x), d);
3349 : }
3350 :
3351 192251 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3352 14096 : case t_FFELT: return FF_inv(x);
3353 5323752 : case t_POL: return gred_rfrac_simple(gen_1,x);
3354 30058 : case t_SER: return ser_inv(x);
3355 2933 : case t_RFRAC:
3356 : {
3357 2933 : GEN n = gel(x,1), d = gel(x,2);
3358 2933 : pari_sp av = avma, ltop;
3359 2933 : if (gequal0(n)) pari_err_INV("ginv",x);
3360 :
3361 2933 : n = simplify_shallow(n);
3362 2933 : if (typ(n) != t_POL || varn(n) != varn(d))
3363 : {
3364 2933 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3365 679 : ltop = avma;
3366 679 : z = RgX_Rg_div(d,n);
3367 : } else {
3368 0 : ltop = avma;
3369 0 : z = cgetg(3,t_RFRAC);
3370 0 : gel(z,1) = RgX_copy(d);
3371 0 : gel(z,2) = RgX_copy(n);
3372 : }
3373 679 : stackdummy(av, ltop);
3374 679 : return z;
3375 : }
3376 :
3377 21 : case t_VEC: if (!is_ext_qfr(x)) break;
3378 : case t_QFB:
3379 153955 : return qfbpow(x, gen_m1);
3380 39581 : case t_MAT:
3381 39581 : y = RgM_inv(x);
3382 39574 : if (!y) pari_err_INV("ginv",x);
3383 39504 : return y;
3384 28 : case t_VECSMALL:
3385 : {
3386 28 : long i, lx = lg(x)-1;
3387 28 : y = zero_zv(lx);
3388 112 : for (i=1; i<=lx; i++)
3389 : {
3390 84 : long xi = x[i];
3391 84 : if (xi<1 || xi>lx || y[xi])
3392 0 : pari_err_TYPE("ginv [not a permutation]", x);
3393 84 : y[xi] = i;
3394 : }
3395 28 : return y;
3396 : }
3397 : }
3398 6 : pari_err_TYPE("inverse",x);
3399 : return NULL; /* LCOV_EXCL_LINE */
3400 : }
|
