Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** GENERIC OPERATIONS **/
18 : /** (first part) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_mod
25 :
26 : /* assume z[1] was created last */
27 : #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
28 : z = gerepileupto((pari_sp)(z+3), gel(z,1));
29 :
30 : /* assume z[1] was created last */
31 : #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
32 : z = gerepileupto((pari_sp)(z+3), gel(z,1));\
33 : else\
34 : gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
35 :
36 : static void
37 105 : warn_coercion(GEN x, GEN y, GEN z)
38 : {
39 105 : if (DEBUGLEVEL)
40 56 : pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
41 105 : }
42 :
43 : static long
44 28 : kro_quad(GEN x, GEN y)
45 28 : { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
46 :
47 : /* is -1 not a square in Zp, assume p prime */
48 : INLINE int
49 42 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
50 :
51 : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
52 : static GEN mulpp(GEN x, GEN y);
53 : static GEN divpp(GEN x, GEN y);
54 : /* Argument codes for inline routines
55 : * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
56 : * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
57 : * T: some type (to be converted to PADIC)
58 : */
59 : static GEN
60 303857536 : addRc(GEN x, GEN y) {
61 303857536 : GEN z = cgetg(3,t_COMPLEX);
62 303842290 : gel(z,1) = gadd(x,gel(y,1));
63 303858117 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 351221376 : mulRc(GEN x, GEN y) {
67 351221376 : GEN z = cgetg(3,t_COMPLEX);
68 351154055 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
69 351408103 : gel(z,2) = gmul(x,gel(y,2)); return z;
70 : }
71 : /* for INTMODs: can't simplify when Re(x) = gen_0 */
72 : static GEN
73 49 : mulRc_direct(GEN x, GEN y) {
74 49 : GEN z = cgetg(3,t_COMPLEX);
75 49 : gel(z,1) = gmul(x,gel(y,1));
76 49 : gel(z,2) = gmul(x,gel(y,2)); return z;
77 : }
78 : static GEN
79 2352679 : divRc(GEN x, GEN y) {
80 2352679 : GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
81 2352670 : GEN z = cgetg(3,t_COMPLEX);
82 2352668 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
83 2352657 : gel(z,2) = gmul(mt, gel(y,2));
84 2352663 : return z;
85 : }
86 : static GEN
87 23109338 : divcR(GEN x, GEN y) {
88 23109338 : GEN z = cgetg(3,t_COMPLEX);
89 23109192 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
90 23109014 : gel(z,2) = gdiv(gel(x,2), y); return z;
91 : }
92 : static GEN
93 1274 : addRq(GEN x, GEN y) {
94 1274 : GEN z = cgetg(4,t_QUAD);
95 1274 : gel(z,1) = ZX_copy(gel(y,1));
96 1274 : gel(z,2) = gadd(x, gel(y,2));
97 1274 : gel(z,3) = gcopy(gel(y,3)); return z;
98 : }
99 : static GEN
100 2352 : mulRq(GEN x, GEN y) {
101 2352 : GEN z = cgetg(4,t_QUAD);
102 2352 : gel(z,1) = ZX_copy(gel(y,1));
103 2352 : gel(z,2) = gmul(x,gel(y,2));
104 2352 : gel(z,3) = gmul(x,gel(y,3)); return z;
105 : }
106 : static GEN
107 77 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
108 77 : long i = gexpo(x) - gexpo(y);
109 77 : if (i > 0) prec += nbits2extraprec(i);
110 77 : return gerepileupto(av, gadd(y, quadtofp(x, prec)));
111 : }
112 : static GEN
113 37507586 : mulrfrac(GEN x, GEN y)
114 : {
115 : pari_sp av;
116 37507586 : GEN z, a = gel(y,1), b = gel(y,2);
117 37507586 : if (is_pm1(a)) /* frequent special case */
118 : {
119 12651453 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
120 12652955 : return z;
121 : }
122 24856062 : av = avma;
123 24856062 : return gerepileuptoleaf(av, divri(mulri(x,a), b));
124 : }
125 : static GEN
126 28 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
127 28 : return gerepileupto(av, gmul(y, quadtofp(x, prec)));
128 : }
129 : static GEN
130 63 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
131 63 : return gerepileupto(av, gdiv(quadtofp(x,prec), y));
132 : }
133 : static GEN
134 42 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
135 42 : return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
136 : }
137 : /* y PADIC, x + y by converting x to padic */
138 : static GEN
139 7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
140 :
141 7 : if (!valp(y)) z = cvtop2(x,y);
142 : else {
143 7 : long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
144 7 : z = cvtop(x, gel(y,2), l);
145 : }
146 7 : return gerepileupto(av, addsub_pp(z, y, addii));
147 : }
148 : /* y PADIC, x * y by converting x to padic */
149 : static GEN
150 4982911 : mulTp(GEN x, GEN y) { pari_sp av = avma;
151 4982911 : return gerepileupto(av, mulpp(cvtop2(x,y), y));
152 : }
153 : /* y PADIC, non zero x / y by converting x to padic */
154 : static GEN
155 3660 : divTp(GEN x, GEN y) { pari_sp av = avma;
156 3660 : return gerepileupto(av, divpp(cvtop2(x,y), y));
157 : }
158 : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
159 : * converted to O(p^e) and division by 0 */
160 : static GEN
161 1204910 : divpT(GEN x, GEN y) { pari_sp av = avma;
162 1204910 : return gerepileupto(av, divpp(x, cvtop2(y,x)));
163 : }
164 :
165 : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
166 : * clean memory from z on */
167 : static GEN
168 3690733 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
169 3690733 : if (lgefint(X) == 3) {
170 3367566 : ulong u = Fl_add(itou(x),itou(y), X[2]);
171 3367566 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
172 : }
173 : else {
174 323167 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
175 323167 : gel(z,2) = gerepileuptoint((pari_sp)z, u);
176 : }
177 3690733 : gel(z,1) = icopy(X); return z;
178 : }
179 : static GEN
180 1158553 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 1158553 : if (lgefint(X) == 3) {
182 784519 : ulong u = Fl_sub(itou(x),itou(y), X[2]);
183 784519 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
184 : }
185 : else {
186 374034 : GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
187 374034 : gel(z,2) = gerepileuptoint((pari_sp)z, u);
188 : }
189 1158553 : gel(z,1) = icopy(X); return z;
190 : }
191 : /* cf add_intmod_same */
192 : static GEN
193 3005196 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
194 3005196 : if (lgefint(X) == 3) {
195 2863550 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
196 2863550 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
197 : }
198 : else
199 141646 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
200 3005197 : gel(z,1) = icopy(X); return z;
201 : }
202 : /* cf add_intmod_same */
203 : static GEN
204 341964 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
205 : {
206 341964 : if (lgefint(X) == 3) {
207 319422 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
208 319359 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 22542 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
212 341901 : gel(z,1) = icopy(X); return z;
213 : }
214 :
215 : /*******************************************************************/
216 : /* */
217 : /* REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC) */
218 : /* */
219 : /* (static routines are not memory clean, but OK for gerepileupto) */
220 : /*******************************************************************/
221 : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
222 : * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
223 : static GEN
224 9857237 : rfrac_denom_mul_scal(GEN d, GEN y)
225 : {
226 9857237 : GEN D = RgX_Rg_mul(d, y);
227 9857237 : if (lg(D) != lg(d))
228 : { /* try to generate a meaningful diagnostic */
229 0 : D = gdiv(leading_coeff(d), y); /* should fail */
230 0 : pari_err_INV("gred_rfrac", y); /* better than nothing */
231 : }
232 9857237 : return D;
233 : }
234 :
235 : static int
236 57829141 : Leading_is_neg(GEN x)
237 : {
238 122278619 : while (typ(x) == t_POL) x = leading_coeff(x);
239 57829141 : return is_real_t(typ(x))? gsigne(x) < 0: 0;
240 : }
241 :
242 : static int
243 154482632 : transtype(GEN x) { return x != gen_1 && typ(x) != t_PADIC; }
244 :
245 : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
246 : GEN
247 57833391 : gred_rfrac_simple(GEN n, GEN d)
248 : {
249 : GEN _1n, _1d, c, cn, cd, z;
250 57833391 : long dd = degpol(d);
251 :
252 57833391 : if (dd <= 0)
253 : {
254 4250 : if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
255 4250 : n = gdiv(n, gel(d,2));
256 4250 : if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
257 4250 : return n;
258 : }
259 57829141 : if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
260 57829141 : _1n = Rg_get_1(n);
261 57829141 : _1d = Rg_get_1(d);
262 57829141 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
263 57829141 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
264 57829134 : cd = content(d);
265 59699609 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
266 57829134 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
267 57829134 : if (!gequal1(cd)) {
268 6568093 : d = RgX_Rg_div(d,cd);
269 6568093 : if (!gequal1(cn))
270 : {
271 1294238 : if (gequal0(cn)) {
272 49 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
273 0 : n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
274 0 : c = gen_1;
275 : } else {
276 1294189 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
277 1294189 : c = gdiv(cn,cd);
278 : }
279 : }
280 : else
281 5273855 : c = ginv(cd);
282 : } else {
283 51261041 : if (!gequal1(cn))
284 : {
285 3277487 : if (gequal0(cn)) {
286 1484 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
287 21 : c = gen_1;
288 : } else {
289 3276003 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
290 3276003 : c = cn;
291 : }
292 : } else {
293 47983554 : GEN y = cgetg(3,t_RFRAC);
294 47983554 : gel(y,1) = gcopy(n);
295 47983554 : gel(y,2) = RgX_copy(d); return y;
296 : }
297 : }
298 :
299 9844068 : if (typ(c) == t_POL)
300 : {
301 900639 : z = c;
302 939069 : do { z = content(z); } while (typ(z) == t_POL);
303 900639 : cd = denom_i(z);
304 900639 : cn = gmul(c, cd);
305 : }
306 : else
307 : {
308 8943429 : cn = numer_i(c);
309 8943429 : cd = denom_i(c);
310 : }
311 9844068 : z = cgetg(3,t_RFRAC);
312 9844068 : gel(z,1) = gmul(n, cn);
313 9844068 : gel(z,2) = d = rfrac_denom_mul_scal(d, cd);
314 : /* can happen: Pol(O(17^-1)) / Pol([Mod(9,23), O(23^-3)]) */
315 9844068 : if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
316 9844068 : return z;
317 : }
318 :
319 : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
320 : static GEN
321 248374 : fix_rfrac(GEN x, long d)
322 : {
323 : GEN z, N, D;
324 248374 : if (!d || typ(x) == t_POL) return x;
325 171710 : z = cgetg(3, t_RFRAC);
326 171710 : N = gel(x,1);
327 171710 : D = gel(x,2);
328 171710 : if (d > 0) {
329 8722 : gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
330 165935 : : monomialcopy(N,d,varn(D));
331 165872 : gel(z, 2) = RgX_copy(D);
332 : } else {
333 5838 : gel(z, 1) = gcopy(N);
334 5838 : gel(z, 2) = RgX_shift(D, -d);
335 : }
336 171710 : return z;
337 : }
338 :
339 : /* assume d != 0 */
340 : static GEN
341 44766283 : gred_rfrac2(GEN n, GEN d)
342 : {
343 : GEN y, z, _1n, _1d;
344 : long v, vd, vn;
345 :
346 44766283 : n = simplify_shallow(n);
347 44766283 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
348 37937776 : d = simplify_shallow(d);
349 37937776 : if (typ(d) != t_POL) return gdiv(n,d);
350 36758990 : vd = varn(d);
351 36758990 : if (typ(n) != t_POL)
352 : {
353 20623729 : if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
354 20622322 : if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
355 0 : pari_err_BUG("gred_rfrac2 [incompatible variables]");
356 : }
357 16135261 : vn = varn(n);
358 16135261 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
359 15997552 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
360 15837779 : _1n = Rg_get_1(n);
361 15837779 : _1d = Rg_get_1(d);
362 15837779 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
363 15837779 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
364 :
365 : /* now n and d are t_POLs in the same variable */
366 15837779 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
367 15837779 : if (!degpol(d))
368 : {
369 12757958 : n = RgX_Rg_div(n,gel(d,2));
370 12757958 : return v? RgX_mulXn(n,v): n;
371 : }
372 :
373 : /* X does not divide gcd(n,d), deg(d) > 0 */
374 3079821 : if (!isinexact(n) && !isinexact(d))
375 : {
376 3079576 : y = RgX_divrem(n, d, &z);
377 3079576 : if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
378 248129 : z = RgX_gcd(d, z);
379 248129 : if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
380 : }
381 248374 : return fix_rfrac(gred_rfrac_simple(n,d), v);
382 : }
383 :
384 : /* x,y t_INT, return x/y in reduced form */
385 : GEN
386 134839700 : Qdivii(GEN x, GEN y)
387 : {
388 134839700 : pari_sp av = avma;
389 : GEN r, q;
390 :
391 134839700 : if (lgefint(y) == 3)
392 : {
393 116991813 : q = Qdiviu(x, y[2]);
394 116988365 : if (signe(y) > 0) return q;
395 11012033 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
396 11012246 : return q;
397 : }
398 17847887 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
399 17849296 : if (equali1(x))
400 : {
401 5178187 : if (!signe(y)) pari_err_INV("gdiv",y);
402 5178082 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
403 : }
404 12671116 : q = dvmdii(x,y,&r);
405 12671008 : if (r == gen_0) return q; /* gen_0 intended */
406 8318273 : r = gcdii(y, r);
407 8318260 : if (lgefint(r) == 3)
408 : {
409 7445183 : ulong t = r[2];
410 7445183 : set_avma(av);
411 7445154 : if (t == 1) q = mkfraccopy(x,y);
412 : else
413 : {
414 2649234 : q = cgetg(3,t_FRAC);
415 2649335 : gel(q,1) = diviuexact(x,t);
416 2649288 : gel(q,2) = diviuexact(y,t);
417 : }
418 : }
419 : else
420 : { /* rare: r and q left on stack for efficiency */
421 873077 : q = cgetg(3,t_FRAC);
422 873083 : gel(q,1) = diviiexact(x,r);
423 873083 : gel(q,2) = diviiexact(y,r);
424 : }
425 8318181 : normalize_frac(q); return q;
426 : }
427 :
428 : /* x t_INT, return x/y in reduced form */
429 : GEN
430 136306783 : Qdiviu(GEN x, ulong y)
431 : {
432 136306783 : pari_sp av = avma;
433 : ulong r, t;
434 : GEN q;
435 :
436 136306783 : if (y == 1) return icopy(x);
437 115520352 : if (!y) pari_err_INV("gdiv",gen_0);
438 115520352 : if (equali1(x)) retmkfrac(gen_1, utoipos(y));
439 108356540 : q = absdiviu_rem(x,y,&r);
440 108352421 : if (!r)
441 : {
442 62511262 : if (signe(x) < 0) togglesign(q);
443 62511446 : return q;
444 : }
445 45841159 : t = ugcd(y, r); set_avma(av);
446 45846963 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
447 18373699 : retmkfrac(diviuexact(x,t), utoipos(y / t));
448 : }
449 :
450 : /* x t_INT, return x/y in reduced form */
451 : GEN
452 1560427 : Qdivis(GEN x, long y)
453 : {
454 1560427 : pari_sp av = avma;
455 : ulong r, t;
456 : long s;
457 : GEN q;
458 :
459 1560427 : if (y > 0) return Qdiviu(x, y);
460 162429 : if (!y) pari_err_INV("gdiv",gen_0);
461 162429 : s = signe(x);
462 162429 : if (!s) return gen_0;
463 132202 : if (y < 0) { y = -y; s = -s; }
464 132202 : if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
465 131922 : if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
466 129759 : q = absdiviu_rem(x,y,&r);
467 129759 : if (!r)
468 : {
469 55356 : if (s < 0) togglesign(q);
470 55356 : return q;
471 : }
472 74403 : t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
473 74403 : if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
474 74403 : gel(q,1) = x; setsigne(x, s);
475 74403 : gel(q,2) = utoipos(y); return q;
476 : }
477 :
478 : /*******************************************************************/
479 : /* */
480 : /* CONJUGATION */
481 : /* */
482 : /*******************************************************************/
483 : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
484 : static GEN
485 18326 : quad_polmod_conj(GEN x, GEN y)
486 : {
487 : GEN z, u, v, a, b;
488 18326 : if (typ(x) != t_POL) return gcopy(x);
489 18326 : if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
490 18326 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
491 18326 : b = gel(y,3); v = gel(x,2);
492 18326 : z = cgetg(4, t_POL); z[1] = x[1];
493 18326 : gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
494 18326 : gel(z,3) = gneg(u); return z;
495 : }
496 : static GEN
497 18326 : quad_polmod_norm(GEN x, GEN y)
498 : {
499 : GEN z, u, v, a, b, c;
500 18326 : if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
501 18326 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
502 18326 : b = gel(y,3); v = gel(x,2);
503 18326 : c = gel(y,2);
504 18326 : z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
505 18326 : if (!gequal1(a)) z = gdiv(z, a);
506 18326 : return gadd(z, gsqr(v));
507 : }
508 :
509 : GEN
510 19573237 : conj_i(GEN x)
511 : {
512 19573237 : switch(typ(x))
513 : {
514 6536612 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
515 6536612 : return x;
516 :
517 12873013 : case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
518 :
519 924 : case t_QUAD:
520 : {
521 924 : GEN y = cgetg(4,t_QUAD);
522 924 : gel(y,1) = gel(x,1);
523 924 : gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
524 924 : : gadd(gel(x,2), gel(x,3));
525 924 : gel(y,3) = gneg(gel(x,3)); return y;
526 : }
527 350 : case t_POL: pari_APPLY_pol_normalized(conj_i(gel(x,i)));
528 31591 : case t_SER: pari_APPLY_ser_normalized(conj_i(gel(x,i)));
529 :
530 152782 : case t_RFRAC:
531 : case t_VEC:
532 : case t_COL:
533 550463 : case t_MAT: pari_APPLY_same(conj_i(gel(x,i)));
534 :
535 0 : case t_POLMOD:
536 : {
537 0 : GEN X = gel(x,1);
538 0 : long d = degpol(X);
539 0 : if (d < 2) return x;
540 0 : if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
541 : }
542 : }
543 0 : pari_err_TYPE("gconj",x);
544 : return NULL; /* LCOV_EXCL_LINE */
545 : }
546 : GEN
547 133181 : gconj(GEN x)
548 133181 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
549 :
550 : GEN
551 84 : conjvec(GEN x,long prec)
552 : {
553 : long lx, s, i;
554 : GEN z;
555 :
556 84 : switch(typ(x))
557 : {
558 0 : case t_INT: case t_INTMOD: case t_FRAC:
559 0 : return mkcolcopy(x);
560 :
561 0 : case t_COMPLEX: case t_QUAD:
562 0 : z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
563 :
564 28 : case t_FFELT:
565 28 : return FF_conjvec(x);
566 :
567 0 : case t_VEC: case t_COL:
568 0 : lx = lg(x); z = cgetg(lx,t_MAT);
569 0 : if (lx == 1) return z;
570 0 : gel(z,1) = conjvec(gel(x,1),prec);
571 0 : s = lgcols(z);
572 0 : for (i=2; i<lx; i++)
573 : {
574 0 : gel(z,i) = conjvec(gel(x,i),prec);
575 0 : if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
576 : }
577 0 : break;
578 :
579 56 : case t_POLMOD: {
580 56 : GEN T = gel(x,1), r;
581 : pari_sp av;
582 :
583 56 : lx = lg(T);
584 56 : if (lx <= 3) return cgetg(1,t_COL);
585 56 : x = gel(x,2);
586 238 : for (i=2; i<lx; i++)
587 : {
588 189 : GEN c = gel(T,i);
589 189 : switch(typ(c)) {
590 7 : case t_INTMOD: {
591 7 : GEN p = gel(c,1);
592 : pari_sp av;
593 7 : if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
594 7 : av = avma;
595 7 : T = RgX_to_FpX(T,p);
596 7 : x = RgX_to_FpX(x, p);
597 7 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
598 7 : z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
599 7 : return gerepileupto(av, z);
600 : }
601 182 : case t_INT:
602 182 : case t_FRAC: break;
603 0 : default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
604 : }
605 : }
606 49 : if (typ(x) != t_POL)
607 : {
608 0 : if (!is_rational_t(typ(x)))
609 0 : pari_err_TYPE("conjvec [not a rational t_POL]",x);
610 0 : retconst_col(lx-3, gcopy(x));
611 : }
612 49 : RgX_check_QX(x,"conjvec");
613 49 : av = avma;
614 49 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
615 49 : r = cleanroots(T,prec);
616 49 : z = cgetg(lx-2,t_COL);
617 182 : for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
618 49 : return gerepileupto(av, z);
619 : }
620 :
621 0 : default:
622 0 : pari_err_TYPE("conjvec",x);
623 : return NULL; /* LCOV_EXCL_LINE */
624 : }
625 0 : return z;
626 : }
627 :
628 : /********************************************************************/
629 : /** **/
630 : /** ADDITION **/
631 : /** **/
632 : /********************************************************************/
633 : /* x, y compatible PADIC, op = add or sub */
634 : static GEN
635 19968435 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
636 : {
637 19968435 : pari_sp av = avma;
638 : long d,e,r,rx,ry;
639 19968435 : GEN u, z, mod, p = gel(x,2);
640 : int swap;
641 :
642 19968435 : (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
643 19967987 : e = valp(x);
644 19967987 : r = valp(y); d = r-e;
645 19967987 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
646 19967987 : rx = precp(x);
647 19967987 : ry = precp(y);
648 19967987 : if (d) /* v(x) < v(y) */
649 : {
650 10703444 : r = d+ry; z = powiu(p,d);
651 10703570 : if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
652 10703592 : z = mulii(z,gel(y,4));
653 10703536 : u = swap? op(z, gel(x,4)): op(gel(x,4), z);
654 : }
655 : else
656 : {
657 : long c;
658 9264543 : if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
659 9264543 : u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
660 9265697 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
661 : {
662 138304 : set_avma(av); return zeropadic(p, e+r);
663 : }
664 9127098 : if (c)
665 : {
666 3351428 : mod = diviiexact(mod, powiu(p,c));
667 3351429 : r -= c;
668 3351429 : e += c;
669 : }
670 : }
671 19830638 : u = modii(u, mod);
672 19829376 : set_avma(av); z = cgetg(5,t_PADIC);
673 19828609 : z[1] = evalprecp(r) | evalvalp(e);
674 19828261 : gel(z,2) = icopy(p);
675 19828064 : gel(z,3) = icopy(mod);
676 19827923 : gel(z,4) = icopy(u); return z;
677 : }
678 : /* Rg_to_Fp(t_FRAC) without GC */
679 : static GEN
680 226838 : Q_to_Fp(GEN x, GEN p)
681 226838 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
682 : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
683 : static GEN
684 4762818 : addQp(GEN x, GEN y)
685 : {
686 4762818 : pari_sp av = avma;
687 4762818 : long d, r, e, vy = valp(y), py = precp(y);
688 4762818 : GEN z, q, mod, u, p = gel(y,2);
689 :
690 4762818 : e = Q_pvalrem(x, p, &x);
691 4762835 : d = vy - e; r = d + py;
692 4762835 : if (r <= 0) { set_avma(av); return gcopy(y); }
693 4760875 : mod = gel(y,3);
694 4760875 : u = gel(y,4);
695 4760875 : (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
696 :
697 4760857 : if (d > 0)
698 : {
699 1376474 : q = powiu(p,d);
700 1376528 : mod = mulii(mod, q);
701 1376533 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
702 1376533 : u = addii(x, mulii(u, q));
703 : }
704 3384383 : else if (d < 0)
705 : {
706 405332 : q = powiu(p,-d);
707 405332 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
708 405332 : u = addii(u, mulii(x, q));
709 405332 : r = py; e = vy;
710 : }
711 : else
712 : {
713 : long c;
714 2979051 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
715 2979051 : u = addii(u, x);
716 2979057 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
717 : {
718 1204 : set_avma(av); return zeropadic(p,e+r);
719 : }
720 2977848 : if (c)
721 : {
722 970645 : mod = diviiexact(mod, powiu(p,c));
723 970645 : r -= c;
724 970645 : e += c;
725 : }
726 : }
727 4759713 : u = modii(u, mod); set_avma(av);
728 4759688 : z = cgetg(5,t_PADIC);
729 4759627 : z[1] = evalprecp(r) | evalvalp(e);
730 4759615 : gel(z,2) = icopy(p);
731 4759591 : gel(z,3) = icopy(mod);
732 4759563 : gel(z,4) = icopy(u); return z;
733 : }
734 :
735 : /* Mod(x,X) + Mod(y,X) */
736 : #define addsub_polmod_same addsub_polmod_scal
737 : /* Mod(x,X) +/- Mod(y,Y) */
738 : static GEN
739 7203 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
740 : {
741 7203 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
742 7203 : GEN z = cgetg(3,t_POLMOD);
743 7203 : long vx = varn(X), vy = varn(Y);
744 7203 : if (vx==vy) {
745 : pari_sp av;
746 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
747 14 : warn_coercion(X,Y,gel(z,1));
748 14 : gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
749 : }
750 7189 : if (varncmp(vx, vy) < 0)
751 7189 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
752 : else
753 0 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
754 7189 : gel(z,2) = op(x, y); return z;
755 : }
756 : /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
757 : static GEN
758 13426300 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
759 13426300 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
760 :
761 : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
762 : static GEN
763 408094 : add_ser_scal(GEN y, GEN x)
764 : {
765 : long i, v, ly, vy;
766 : GEN z;
767 :
768 408094 : if (isrationalzero(x)) return gcopy(y);
769 381837 : ly = lg(y);
770 381837 : v = valser(y);
771 381837 : if (v < 3-ly) return gcopy(y);
772 : /* v + ly >= 3 */
773 381585 : if (v < 0)
774 : {
775 1162 : z = cgetg(ly,t_SER); z[1] = y[1];
776 3276 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
777 1162 : gel(z,i) = gadd(x,gel(y,i)); i++;
778 3157 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
779 1162 : return normalizeser(z);
780 : }
781 380423 : vy = varn(y);
782 380423 : if (v > 0)
783 : {
784 19299 : if (ser_isexactzero(y))
785 9338 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
786 9961 : y -= v; ly += v;
787 9961 : z = cgetg(ly,t_SER);
788 9961 : x = gcopy(x);
789 20559 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
790 : }
791 361124 : else if (ser_isexactzero(y)) return gcopy(y);
792 : else
793 : { /* v = 0, ly >= 3 */
794 361117 : z = cgetg(ly,t_SER);
795 361117 : x = gadd(x, gel(y,2));
796 361117 : i = 3;
797 : }
798 1598979 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
799 371078 : gel(z,2) = x;
800 371078 : z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
801 371078 : return gequal0(x)? normalizeser(z): z;
802 : }
803 : static long
804 7161978 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
805 : /* x,y t_SER in the same variable: x+y */
806 : static GEN
807 3581388 : ser_add(GEN x, GEN y)
808 : {
809 3581388 : long i, lx,ly, n = valser(y) - valser(x);
810 : GEN z;
811 3581388 : if (n < 0) { n = -n; swap(x,y); }
812 : /* valser(x) <= valser(y) */
813 3581388 : lx = _serprec(x);
814 3581388 : if (lx == 2) /* don't lose type information */
815 : {
816 798 : z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
817 798 : setvalser(z, valser(x)); return z;
818 : }
819 3580590 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
820 3580590 : if (n)
821 : {
822 107782 : if (n+2 > lx) return gcopy(x);
823 107082 : z = cgetg(ly,t_SER);
824 819116 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
825 506624 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
826 : } else {
827 3472808 : z = cgetg(ly,t_SER);
828 17589354 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
829 : }
830 3579890 : z[1] = x[1]; return normalizeser(z);
831 : }
832 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
833 : static GEN
834 8830966 : add_rfrac_scal(GEN y, GEN x)
835 : {
836 : pari_sp av;
837 : GEN n;
838 :
839 8830966 : if (isintzero(x)) return gcopy(y); /* frequent special case */
840 5157194 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
841 5157194 : return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
842 : }
843 :
844 : /* x "scalar", ty != t_MAT and nonscalar */
845 : static GEN
846 26106544 : add_scal(GEN y, GEN x, long ty)
847 : {
848 26106544 : switch(ty)
849 : {
850 21194327 : case t_POL: return RgX_Rg_add(y, x);
851 408066 : case t_SER: return add_ser_scal(y, x);
852 4468785 : case t_RFRAC: return add_rfrac_scal(y, x);
853 0 : case t_COL: return RgC_Rg_add(y, x);
854 35356 : case t_VEC:
855 35356 : if (isintzero(x)) return gcopy(y);
856 175 : break;
857 : }
858 185 : pari_err_TYPE2("+",x,y);
859 : return NULL; /* LCOV_EXCL_LINE */
860 : }
861 :
862 : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
863 : static GEN
864 14927260 : setfrac(GEN z, GEN a, GEN b)
865 : {
866 14927260 : gel(z,1) = icopy_avma(a, (pari_sp)z);
867 14927253 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
868 14927325 : set_avma((pari_sp)gel(z,2)); return z;
869 : }
870 : /* z <- a / (b*Q), (Q,a) = 1 */
871 : static GEN
872 14120677 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
873 : {
874 14120677 : GEN q = Qdivii(a, b); /* != 0 */
875 14120723 : if (typ(q) == t_INT)
876 : {
877 1976933 : gel(z,1) = gerepileuptoint((pari_sp)Q, q);
878 1976933 : gel(z,2) = Q; return z;
879 : }
880 12143790 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
881 : }
882 : static GEN
883 28647506 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
884 : {
885 28647506 : GEN x1 = gel(x,1), x2 = gel(x,2);
886 28647506 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
887 28647506 : int s = cmpii(x2, y2);
888 :
889 : /* common denominator: (x1 op y1) / x2 */
890 28647501 : if (!s)
891 : {
892 10186725 : pari_sp av = avma;
893 10186725 : return gerepileupto(av, Qdivii(op(x1, y1), x2));
894 : }
895 18460776 : z = cgetg(3, t_FRAC);
896 18460789 : if (s < 0)
897 : {
898 10461721 : Q = dvmdii(y2, x2, &r);
899 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
900 10461681 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
901 2762866 : delta = gcdii(x2,r);
902 : }
903 : else
904 : {
905 7999068 : Q = dvmdii(x2, y2, &r);
906 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
907 7999082 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
908 1577179 : delta = gcdii(y2,r);
909 : }
910 : /* delta = gcd(x2,y2) */
911 4340084 : if (equali1(delta))
912 : { /* numerator is nonzero */
913 1556547 : gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
914 1556547 : gel(z,2) = mulii(x2,y2); return z;
915 : }
916 2783533 : x2 = diviiexact(x2,delta);
917 2783537 : y2 = diviiexact(y2,delta);
918 2783535 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
919 2783543 : q = dvmdii(n, delta, &r);
920 2783537 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
921 2528777 : r = gcdii(delta, r);
922 2528783 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
923 2528783 : return setfrac(z, n, mulii(mulii(x2, y2), delta));
924 : }
925 :
926 : /* assume x2, y2 are t_POLs in the same variable */
927 : static GEN
928 3039747 : add_rfrac(GEN x, GEN y)
929 : {
930 3039747 : pari_sp av = avma;
931 3039747 : GEN x1 = gel(x,1), x2 = gel(x,2);
932 3039747 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
933 :
934 3039747 : delta = RgX_gcd(x2,y2);
935 3039747 : if (!degpol(delta))
936 : {
937 672 : n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
938 672 : d = RgX_mul(x2, y2);
939 672 : return gerepileupto(av, gred_rfrac_simple(n, d));
940 : }
941 3039075 : x2 = RgX_div(x2,delta);
942 3039075 : y2 = RgX_div(y2,delta);
943 3039075 : n = gadd(gmul(x1,y2), gmul(y1,x2));
944 3039075 : if (!signe(n))
945 : {
946 721425 : n = simplify_shallow(n);
947 721425 : if (isexactzero(n))
948 : {
949 721418 : if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
950 14 : return gerepileupto(av, scalarpol(n, varn(x2)));
951 : }
952 7 : return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
953 : }
954 2317650 : if (degpol(n) == 0)
955 1150603 : return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
956 1167047 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
957 1167047 : if (isexactzero(r))
958 : {
959 : GEN z;
960 228087 : d = RgX_mul(x2, y2);
961 : /* "constant" denominator ? */
962 228087 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
963 228087 : return gerepileupto(av, z);
964 : }
965 938960 : r = RgX_gcd(delta, r);
966 938960 : if (degpol(r))
967 : {
968 160740 : n = RgX_div(n, r);
969 160740 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
970 : }
971 : else
972 778220 : d = RgX_mul(gel(x,2), y2);
973 938960 : return gerepileupto(av, gred_rfrac_simple(n, d));
974 : }
975 :
976 : GEN
977 5883807095 : gadd(GEN x, GEN y)
978 : {
979 5883807095 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
980 : pari_sp av, tetpil;
981 : GEN z, p1;
982 :
983 5883807095 : if (tx == ty) switch(tx) /* shortcut to generic case */
984 : {
985 2911611359 : case t_INT: return addii(x,y);
986 1636653017 : case t_REAL: return addrr(x,y);
987 1521962 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
988 1521962 : z = cgetg(3,t_INTMOD);
989 1521962 : if (X==Y || equalii(X,Y))
990 1521948 : return add_intmod_same(z, X, gel(x,2), gel(y,2));
991 14 : gel(z,1) = gcdii(X,Y);
992 14 : warn_coercion(X,Y,gel(z,1));
993 14 : av = avma; p1 = addii(gel(x,2),gel(y,2));
994 14 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
995 : }
996 23196862 : case t_FRAC: return addsub_frac(x,y,addii);
997 354328617 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
998 353840003 : gel(z,2) = gadd(gel(x,2),gel(y,2));
999 354851006 : if (isintzero(gel(z,2)))
1000 : {
1001 521178 : set_avma((pari_sp)(z+3));
1002 521178 : return gadd(gel(x,1),gel(y,1));
1003 : }
1004 354155137 : gel(z,1) = gadd(gel(x,1),gel(y,1));
1005 354280093 : return z;
1006 14530828 : case t_PADIC:
1007 14530828 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
1008 14530640 : return addsub_pp(x,y, addii);
1009 672 : case t_QUAD: z = cgetg(4,t_QUAD);
1010 672 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1011 672 : gel(z,1) = ZX_copy(gel(x,1));
1012 672 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1013 672 : gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
1014 10161232 : case t_POLMOD:
1015 10161232 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1016 10154063 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
1017 7168 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
1018 14006809 : case t_FFELT: return FF_add(x,y);
1019 72966173 : case t_POL:
1020 72966173 : vx = varn(x);
1021 72966173 : vy = varn(y);
1022 72966173 : if (vx != vy) {
1023 855113 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1024 66302 : else return RgX_Rg_add(y, x);
1025 : }
1026 72111060 : return RgX_add(x, y);
1027 3576033 : case t_SER:
1028 3576033 : vx = varn(x);
1029 3576033 : vy = varn(y);
1030 3576033 : if (vx != vy) {
1031 28 : if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1032 21 : else return add_ser_scal(y, x);
1033 : }
1034 3576005 : return ser_add(x, y);
1035 4331648 : case t_RFRAC:
1036 4331648 : vx = varn(gel(x,2));
1037 4331648 : vy = varn(gel(y,2));
1038 4331648 : if (vx != vy) {
1039 1291901 : if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1040 538397 : else return add_rfrac_scal(y, x);
1041 : }
1042 3039747 : return add_rfrac(x,y);
1043 2060573 : case t_VEC:
1044 2060573 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1045 2060573 : return RgV_add(x,y);
1046 1108619 : case t_COL:
1047 1108619 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1048 1108619 : return RgC_add(x,y);
1049 671356 : case t_MAT:
1050 671356 : lx = lg(x);
1051 671356 : if (lg(y) != lx) pari_err_OP("+",x,y);
1052 671356 : if (lx == 1) return cgetg(1, t_MAT);
1053 671356 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1054 671349 : return RgM_add(x,y);
1055 :
1056 0 : default: pari_err_TYPE2("+",x,y);
1057 : }
1058 : /* tx != ty */
1059 837451782 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1060 :
1061 837451782 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1062 : {
1063 740852540 : case t_INT:
1064 : switch(ty)
1065 : {
1066 436538655 : case t_REAL: return addir(x,y);
1067 2139209 : case t_INTMOD:
1068 2139209 : z = cgetg(3, t_INTMOD);
1069 2139209 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1070 45240460 : case t_FRAC: z = cgetg(3,t_FRAC);
1071 45240464 : gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1072 45240407 : gel(z,2) = icopy(gel(y,2)); return z;
1073 250031965 : case t_COMPLEX: return addRc(x, y);
1074 5596802 : case t_PADIC:
1075 5596802 : if (!signe(x)) return gcopy(y);
1076 4533077 : return addQp(x,y);
1077 1113 : case t_QUAD: return addRq(x, y);
1078 1357806 : case t_FFELT: return FF_Z_add(y,x);
1079 : }
1080 :
1081 : case t_REAL:
1082 57071374 : switch(ty)
1083 : {
1084 13521329 : case t_FRAC:
1085 13521329 : if (!signe(gel(y,1))) return rcopy(x);
1086 13521329 : if (!signe(x))
1087 : {
1088 12089 : lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1089 12089 : return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1090 : }
1091 13509240 : av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
1092 13508326 : return gerepile(av,tetpil,divri(z,gel(y,2)));
1093 43549975 : case t_COMPLEX: return addRc(x, y);
1094 70 : case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, realprec(x));
1095 :
1096 0 : default: pari_err_TYPE2("+",x,y);
1097 : }
1098 :
1099 17640 : case t_INTMOD:
1100 : switch(ty)
1101 : {
1102 17507 : case t_FRAC: { GEN X = gel(x,1);
1103 17507 : z = cgetg(3, t_INTMOD);
1104 17507 : p1 = Fp_div(gel(y,1), gel(y,2), X);
1105 17507 : return add_intmod_same(z, X, p1, gel(x,2));
1106 : }
1107 14 : case t_FFELT:
1108 14 : if (!equalii(gel(x,1),FF_p_i(y)))
1109 0 : pari_err_OP("+",x,y);
1110 14 : return FF_Z_add(y,gel(x,2));
1111 91 : case t_COMPLEX: return addRc(x, y);
1112 0 : case t_PADIC: { GEN X = gel(x,1);
1113 0 : z = cgetg(3, t_INTMOD);
1114 0 : return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1115 : }
1116 28 : case t_QUAD: return addRq(x, y);
1117 : }
1118 :
1119 : case t_FRAC:
1120 : switch (ty)
1121 : {
1122 10299312 : case t_COMPLEX: return addRc(x, y);
1123 229736 : case t_PADIC:
1124 229736 : if (!signe(gel(x,1))) return gcopy(y);
1125 229736 : return addQp(x,y);
1126 133 : case t_QUAD: return addRq(x, y);
1127 1337 : case t_FFELT: return FF_Q_add(y, x);
1128 : }
1129 :
1130 : case t_FFELT:
1131 0 : pari_err_TYPE2("+",x,y);
1132 :
1133 35 : case t_COMPLEX:
1134 : switch(ty)
1135 : {
1136 28 : case t_PADIC:
1137 28 : return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
1138 7 : case t_QUAD:
1139 7 : lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1140 7 : return gequal0(y)? gcopy(x): addqf(y, x, lx);
1141 : }
1142 :
1143 : case t_PADIC: /* ty == t_QUAD */
1144 0 : return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
1145 : }
1146 : /* tx < ty, !is_const_t(y) */
1147 31818859 : switch(ty)
1148 : {
1149 7951 : case t_MAT:
1150 7951 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1151 7951 : if (isrationalzero(x)) return gcopy(y);
1152 7874 : return RgM_Rg_add(y, x);
1153 206994 : case t_COL:
1154 206994 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1155 206994 : return RgC_Rg_add(y, x);
1156 2393427 : case t_POLMOD: /* is_const_t(tx) in this case */
1157 2393427 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1158 : }
1159 29210487 : if (is_scalar_t(tx)) {
1160 26111554 : if (tx == t_POLMOD)
1161 : {
1162 112774 : vx = varn(gel(x,1));
1163 112774 : vy = gvar(y);
1164 112774 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1165 : else
1166 85544 : if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1167 27755 : return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1168 : }
1169 25998780 : return add_scal(y, x, ty);
1170 : }
1171 : /* x and y are not scalars, ty != t_MAT */
1172 3098940 : vx = gvar(x);
1173 3098947 : vy = gvar(y);
1174 3098947 : if (vx != vy) { /* x or y is treated as a scalar */
1175 22759 : if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1176 32417 : return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1177 32417 : : add_scal(y, x, ty);
1178 : }
1179 : /* vx = vy */
1180 3076188 : switch(tx)
1181 : {
1182 3075691 : case t_POL:
1183 : switch (ty)
1184 : {
1185 5411 : case t_SER:
1186 5411 : if (lg(x) == 2) return gcopy(y);
1187 5390 : i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1188 5390 : i = lg(y) + valser(y) - i;
1189 5390 : if (i < 3) return gcopy(y);
1190 5383 : p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1191 5383 : settyp(p1, t_VECSMALL); /* p1 left on stack */
1192 5383 : return y;
1193 :
1194 3070280 : case t_RFRAC: return add_rfrac_scal(y, x);
1195 : }
1196 0 : break;
1197 :
1198 497 : case t_SER:
1199 497 : if (ty == t_RFRAC)
1200 : {
1201 : long vn, vd;
1202 497 : av = avma;
1203 497 : vn = gval(gel(y,1), vy);
1204 497 : vd = RgX_valrem_inexact(gel(y,2), NULL);
1205 497 : if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
1206 :
1207 497 : l = lg(x) + valser(x) - (vn - vd);
1208 497 : if (l < 3) { set_avma(av); return gcopy(x); }
1209 497 : return gerepileupto(av, gadd(x, rfrac_to_ser_i(y, l)));
1210 : }
1211 0 : break;
1212 : }
1213 0 : pari_err_TYPE2("+",x,y);
1214 : return NULL; /* LCOV_EXCL_LINE */
1215 : }
1216 :
1217 : GEN
1218 270571147 : gaddsg(long x, GEN y)
1219 : {
1220 270571147 : long ty = typ(y);
1221 : GEN z;
1222 :
1223 270571147 : switch(ty)
1224 : {
1225 124841516 : case t_INT: return addsi(x,y);
1226 119175574 : case t_REAL: return addsr(x,y);
1227 12013 : case t_INTMOD:
1228 12013 : z = cgetg(3, t_INTMOD);
1229 12013 : return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1230 15247663 : case t_FRAC: z = cgetg(3,t_FRAC);
1231 15247663 : gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1232 15247663 : gel(z,2) = icopy(gel(y,2)); return z;
1233 8164229 : case t_COMPLEX:
1234 8164229 : z = cgetg(3, t_COMPLEX);
1235 8164229 : gel(z,1) = gaddsg(x, gel(y,1));
1236 8164229 : gel(z,2) = gcopy(gel(y,2)); return z;
1237 :
1238 3130152 : default: return gadd(stoi(x), y);
1239 : }
1240 : }
1241 :
1242 : GEN
1243 3167508 : gsubsg(long x, GEN y)
1244 : {
1245 : GEN z, a, b;
1246 : pari_sp av;
1247 :
1248 3167508 : switch(typ(y))
1249 : {
1250 276674 : case t_INT: return subsi(x,y);
1251 1224993 : case t_REAL: return subsr(x,y);
1252 56 : case t_INTMOD:
1253 56 : z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1254 56 : return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1255 732332 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1256 732332 : gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
1257 732332 : gel(z,2) = icopy(gel(y,2)); return z;
1258 896339 : case t_COMPLEX:
1259 896339 : z = cgetg(3, t_COMPLEX);
1260 896339 : gel(z,1) = gsubsg(x, gel(y,1));
1261 896339 : gel(z,2) = gneg(gel(y,2)); return z;
1262 : }
1263 37114 : av = avma;
1264 37114 : return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
1265 : }
1266 :
1267 : /********************************************************************/
1268 : /** **/
1269 : /** SUBTRACTION **/
1270 : /** **/
1271 : /********************************************************************/
1272 :
1273 : GEN
1274 2828572873 : gsub(GEN x, GEN y)
1275 : {
1276 2828572873 : long tx = typ(x), ty = typ(y);
1277 : pari_sp av;
1278 : GEN z;
1279 2828572873 : if (tx == ty) switch(tx) /* shortcut to generic case */
1280 : {
1281 2059699242 : case t_INT: return subii(x,y);
1282 577204512 : case t_REAL: return subrr(x,y);
1283 1158567 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1284 1158567 : z = cgetg(3,t_INTMOD);
1285 1158567 : if (X==Y || equalii(X,Y))
1286 1158553 : return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1287 14 : gel(z,1) = gcdii(X,Y);
1288 14 : warn_coercion(X,Y,gel(z,1));
1289 14 : av = avma; p1 = subii(gel(x,2),gel(y,2));
1290 14 : gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
1291 : }
1292 5450657 : case t_FRAC: return addsub_frac(x,y, subii);
1293 101975439 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1294 102003486 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1295 101925887 : if (isintzero(gel(z,2)))
1296 : {
1297 21301 : set_avma((pari_sp)(z+3));
1298 21301 : return gsub(gel(x,1),gel(y,1));
1299 : }
1300 101904535 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1301 101893675 : return z;
1302 5437807 : case t_PADIC:
1303 5437807 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
1304 5437807 : return addsub_pp(x,y, subii);
1305 168 : case t_QUAD: z = cgetg(4,t_QUAD);
1306 168 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1307 168 : gel(z,1) = ZX_copy(gel(x,1));
1308 168 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1309 168 : gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1310 851090 : case t_POLMOD:
1311 851090 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1312 851055 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1313 35 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1314 203220 : case t_FFELT: return FF_sub(x,y);
1315 10083788 : case t_POL: {
1316 10083788 : long vx = varn(x);
1317 10083788 : long vy = varn(y);
1318 10083788 : if (vx != vy) {
1319 23622 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1320 2261 : else return Rg_RgX_sub(x, y);
1321 : }
1322 10060166 : return RgX_sub(x, y);
1323 : }
1324 299536 : case t_VEC:
1325 299536 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1326 299536 : return RgV_sub(x,y);
1327 3567164 : case t_COL:
1328 3567164 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1329 3567164 : return RgC_sub(x,y);
1330 251801 : case t_MAT: {
1331 251801 : long lx = lg(x);
1332 251801 : if (lg(y) != lx) pari_err_OP("+",x,y);
1333 251803 : if (lx == 1) return cgetg(1, t_MAT);
1334 172229 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1335 172229 : return RgM_sub(x,y);
1336 : }
1337 2469828 : case t_RFRAC: case t_SER: break;
1338 :
1339 0 : default: pari_err_TYPE2("+",x,y);
1340 : }
1341 64070921 : av = avma;
1342 64070921 : return gerepileupto(av, gadd(x,gneg_i(y)));
1343 : }
1344 :
1345 : /********************************************************************/
1346 : /** **/
1347 : /** MULTIPLICATION **/
1348 : /** **/
1349 : /********************************************************************/
1350 : static GEN
1351 316176 : mul_ser_scal(GEN x, GEN t)
1352 : {
1353 316176 : if (isexactzero(t)) return gmul(Rg_get_0(x), t);
1354 312928 : if (isint1(t)) return gcopy(x);
1355 262360 : if (ser_isexactzero(x))
1356 : {
1357 378 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gmul(gel(x,2), t), varn(x), 1);
1358 378 : setvalser(z, valser(x)); return z;
1359 : }
1360 3442166 : pari_APPLY_ser(gmul(gel(x,i), t));
1361 : }
1362 : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1363 : * [n/d a valid RFRAC] */
1364 : static GEN
1365 10452363 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1366 : {
1367 10452363 : pari_sp av = avma;
1368 : GEN z;
1369 :
1370 10452363 : switch(typ(x))
1371 : {
1372 21 : case t_PADIC:
1373 21 : n = gmul(n, x);
1374 21 : d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
1375 21 : return gerepileupto(av, gdiv(n,d));
1376 :
1377 518 : case t_INTMOD: case t_POLMOD:
1378 518 : n = gmul(n, x);
1379 518 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1380 518 : return gerepileupto(av, gdiv(n,d));
1381 : }
1382 10451824 : z = gred_rfrac2(x, d);
1383 10451824 : n = simplify_shallow(n);
1384 10451824 : if (typ(z) == t_RFRAC)
1385 : {
1386 7936098 : n = gmul(gel(z,1), n);
1387 7936098 : d = gel(z,2);
1388 7936098 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1389 0 : z = RgX_Rg_div(n, d);
1390 : else
1391 7936098 : z = gred_rfrac_simple(n, d);
1392 : }
1393 : else
1394 2515726 : z = gmul(z, n);
1395 10451824 : return gerepileupto(av, z);
1396 : }
1397 : static GEN
1398 122125992 : mul_scal(GEN y, GEN x, long ty)
1399 : {
1400 122125992 : switch(ty)
1401 : {
1402 113232273 : case t_POL:
1403 113232273 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1404 110078801 : return RgX_Rg_mul(y, x);
1405 308035 : case t_SER: return mul_ser_scal(y, x);
1406 8585675 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1407 14 : case t_QFB:
1408 14 : if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1409 : }
1410 12 : pari_err_TYPE2("*",x,y);
1411 : return NULL; /* LCOV_EXCL_LINE */
1412 : }
1413 : static GEN
1414 10542512 : mul_self_scal(GEN x, GEN y)
1415 643179386 : { pari_APPLY_same(gmul(y,gel(x,i))); }
1416 :
1417 : static GEN
1418 160886 : mul_gen_rfrac(GEN X, GEN Y)
1419 : {
1420 160886 : GEN y1 = gel(Y,1), y2 = gel(Y,2);
1421 160886 : long vx = gvar(X), vy = varn(y2);
1422 166647 : return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1423 5761 : gred_rfrac_simple(gmul(y1,X), y2);
1424 : }
1425 : /* (x1/x2) * (y1/y2) */
1426 : static GEN
1427 7907310 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1428 : {
1429 : GEN z, X, Y;
1430 7907310 : pari_sp av = avma;
1431 :
1432 7907310 : X = gred_rfrac2(x1, y2);
1433 7907310 : Y = gred_rfrac2(y1, x2);
1434 7907310 : if (typ(X) == t_RFRAC)
1435 : {
1436 6628116 : if (typ(Y) == t_RFRAC) {
1437 6562519 : x1 = gel(X,1);
1438 6562519 : x2 = gel(X,2);
1439 6562519 : y1 = gel(Y,1);
1440 6562519 : y2 = gel(Y,2);
1441 6562519 : z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1442 : } else
1443 65597 : z = mul_gen_rfrac(Y, X);
1444 : }
1445 1279194 : else if (typ(Y) == t_RFRAC)
1446 95289 : z = mul_gen_rfrac(X, Y);
1447 : else
1448 1183905 : z = gmul(X, Y);
1449 7907310 : return gerepileupto(av, z);
1450 : }
1451 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1452 : static GEN
1453 271003 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1454 : {
1455 271003 : pari_sp av = avma;
1456 271003 : GEN X = gred_rfrac2(x1, y2);
1457 271003 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1458 : {
1459 264423 : x2 = RgX_mul(gel(X,2), x2);
1460 264423 : x1 = gel(X,1);
1461 : }
1462 : else
1463 6580 : x1 = X;
1464 271003 : return gerepileupto(av, gred_rfrac_simple(x1, x2));
1465 : }
1466 :
1467 : /* Mod(y, Y) * x, assuming x scalar */
1468 : static GEN
1469 3427266 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1470 3427266 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1471 :
1472 : /* cf mulqq */
1473 : static GEN
1474 5858691 : quad_polmod_mul(GEN P, GEN x, GEN y)
1475 : {
1476 5858691 : GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
1477 5858691 : pari_sp tetpil, av = avma;
1478 5858691 : T[1] = x[1];
1479 5858691 : p2 = gmul(gel(x,2), gel(y,2));
1480 5858691 : p3 = gmul(gel(x,3), gel(y,3));
1481 5858691 : p1 = gmul(gneg_i(c),p3);
1482 : /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
1483 5858691 : if (typ(b) == t_INT)
1484 : {
1485 5858670 : if (signe(b))
1486 : {
1487 4284101 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1488 4284101 : if (is_pm1(b))
1489 : {
1490 4283443 : if (signe(b) > 0) p3 = gneg(p3);
1491 : }
1492 : else
1493 658 : p3 = gmul(negi(b), p3);
1494 : }
1495 : else
1496 : {
1497 1574569 : p3 = gmul(gel(x,2),gel(y,3));
1498 1574569 : p4 = gmul(gel(x,3),gel(y,2));
1499 : }
1500 : }
1501 : else
1502 : {
1503 21 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1504 21 : p3 = gmul(gneg_i(b), p3);
1505 : }
1506 5858691 : tetpil = avma;
1507 5858691 : gel(T,2) = gadd(p2, p1);
1508 5858691 : gel(T,3) = gadd(p4, p3);
1509 5858691 : gerepilecoeffssp(av,tetpil,T+2,2);
1510 5858691 : return normalizepol_lg(T,4);
1511 : }
1512 : /* Mod(x,T) * Mod(y,T) */
1513 : static GEN
1514 8563693 : mul_polmod_same(GEN T, GEN x, GEN y)
1515 : {
1516 8563693 : GEN z = cgetg(3,t_POLMOD), a;
1517 8563693 : long v = varn(T), lx = lg(x), ly = lg(y);
1518 8563693 : gel(z,1) = RgX_copy(T);
1519 : /* x * y mod T optimised */
1520 8563693 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1521 7892857 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1522 1587375 : a = gmul(x, y);
1523 : else
1524 : {
1525 6976318 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1526 5853420 : a = quad_polmod_mul(T, x, y);
1527 : else
1528 1122898 : a = RgXQ_mul(x, y, gel(z,1));
1529 : }
1530 8563693 : gel(z,2) = a; return z;
1531 : }
1532 : static GEN
1533 484415 : sqr_polmod(GEN T, GEN x)
1534 : {
1535 484415 : GEN a, z = cgetg(3,t_POLMOD);
1536 484415 : gel(z,1) = RgX_copy(T);
1537 484415 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1538 93745 : a = gsqr(x);
1539 : else
1540 : {
1541 390670 : pari_sp av = avma;
1542 390670 : a = RgXQ_sqr(x, gel(z,1));
1543 390670 : a = gerepileupto(av, a);
1544 : }
1545 484415 : gel(z,2) = a; return z;
1546 : }
1547 : /* Mod(x,X) * Mod(y,Y) */
1548 : static GEN
1549 2675106 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1550 : {
1551 2675106 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1552 2675106 : long vx = varn(X), vy = varn(Y);
1553 2675106 : GEN z = cgetg(3,t_POLMOD);
1554 :
1555 2675106 : if (vx==vy) {
1556 : pari_sp av;
1557 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
1558 14 : warn_coercion(X,Y,gel(z,1));
1559 14 : gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
1560 14 : return z;
1561 : }
1562 2675092 : if (varncmp(vx, vy) < 0)
1563 410662 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1564 : else
1565 2264430 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1566 2675092 : gel(z,2) = gmul(x, y); return z;
1567 : }
1568 :
1569 : #if 0 /* used by 3M only */
1570 : /* set z = x+y and return 1 if x,y have the same sign
1571 : * set z = x-y and return 0 otherwise */
1572 : static int
1573 : did_add(GEN x, GEN y, GEN *z)
1574 : {
1575 : long tx = typ(x), ty = typ(y);
1576 : if (tx == ty) switch(tx)
1577 : {
1578 : case t_INT: *z = addii(x,y); return 1;
1579 : case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1580 : case t_REAL:
1581 : if (signe(x) == -signe(y))
1582 : { *z = subrr(x,y); return 0; }
1583 : else
1584 : { *z = addrr(x,y); return 1; }
1585 : }
1586 : if (tx == t_REAL) switch(ty)
1587 : {
1588 : case t_INT:
1589 : if (signe(x) == -signe(y))
1590 : { *z = subri(x,y); return 0; }
1591 : else
1592 : { *z = addri(x,y); return 1; }
1593 : case t_FRAC:
1594 : if (signe(x) == -signe(gel(y,1)))
1595 : { *z = gsub(x,y); return 0; }
1596 : else
1597 : { *z = gadd(x,y); return 1; }
1598 : }
1599 : else if (ty == t_REAL) switch(tx)
1600 : {
1601 : case t_INT:
1602 : if (signe(x) == -signe(y))
1603 : { *z = subir(x,y); return 0; }
1604 : else
1605 : { *z = addir(x,y); return 1; }
1606 : case t_FRAC:
1607 : if (signe(gel(x,1)) == -signe(y))
1608 : { *z = gsub(x,y); return 0; }
1609 : else
1610 : { *z = gadd(x,y); return 1; }
1611 : }
1612 : *z = gadd(x,y); return 1;
1613 : }
1614 : #endif
1615 : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1616 : static GEN
1617 11682598 : mulcIR(GEN x, GEN y)
1618 : {
1619 11682598 : GEN z = cgetg(3,t_COMPLEX);
1620 11682735 : pari_sp av = avma;
1621 11682735 : gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
1622 11682953 : gel(z,2) = gmul(y, gel(x,1));
1623 11682878 : return z;
1624 :
1625 : }
1626 : /* x,y COMPLEX */
1627 : static GEN
1628 273757195 : mulcc(GEN x, GEN y)
1629 : {
1630 273757195 : GEN xr = gel(x,1), xi = gel(x,2);
1631 273757195 : GEN yr = gel(y,1), yi = gel(y,2);
1632 : GEN p1, p2, p3, p4, z;
1633 : pari_sp tetpil, av;
1634 :
1635 273757195 : if (isintzero(xr))
1636 : {
1637 15650107 : if (isintzero(yr)) {
1638 7122769 : av = avma;
1639 7122769 : return gerepileupto(av, gneg(gmul(xi,yi)));
1640 : }
1641 8527155 : return mulcIR(y, xi);
1642 : }
1643 258094319 : if (isintzero(yr)) return mulcIR(x, yi);
1644 :
1645 254932815 : z = cgetg(3,t_COMPLEX); av = avma;
1646 : #if 0
1647 : /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1648 : * e.g. xr + xi if exponents differ */
1649 : if (did_add(xr, xi, &p3))
1650 : {
1651 : if (did_add(yr, yi, &p4)) {
1652 : /* R = xr*yr - xi*yi
1653 : * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1654 : p1 = gmul(xr,yr);
1655 : p2 = gmul(xi,yi); p2 = gneg(p2);
1656 : p3 = gmul(p3, p4);
1657 : p4 = gsub(p2, p1);
1658 : } else {
1659 : /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1660 : * I = xr*yi + xi*yr */
1661 : p1 = gmul(p3,p4);
1662 : p3 = gmul(xr,yi);
1663 : p4 = gmul(xi,yr);
1664 : p2 = gsub(p3, p4);
1665 : }
1666 : } else {
1667 : if (did_add(yr, yi, &p4)) {
1668 : /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1669 : * I = xr*yi +xi*yr */
1670 : p1 = gmul(p3,p4);
1671 : p3 = gmul(xr,yi);
1672 : p4 = gmul(xi,yr);
1673 : p2 = gsub(p4, p3);
1674 : } else {
1675 : /* R = xr*yr - xi*yi
1676 : * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1677 : p3 = gneg( gmul(p3, p4) );
1678 : p1 = gmul(xr,yr);
1679 : p2 = gmul(xi,yi);
1680 : p4 = gadd(p1, p2);
1681 :
1682 : p2 = gneg(p2);
1683 : }
1684 : }
1685 : tetpil = avma;
1686 : gel(z,1) = gadd(p1,p2);
1687 : gel(z,2) = gadd(p3,p4);
1688 : #else
1689 254932571 : if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1690 : { /* 3M formula */
1691 559286 : p3 = addii(xr,xi);
1692 559286 : p4 = addii(yr,yi);
1693 559286 : p1 = mulii(xr,yr);
1694 559286 : p2 = mulii(xi,yi);
1695 559286 : p3 = mulii(p3,p4);
1696 559286 : p4 = addii(p2,p1);
1697 559286 : tetpil = avma;
1698 559286 : gel(z,1) = subii(p1,p2);
1699 559286 : gel(z,2) = subii(p3,p4);
1700 559286 : if (!signe(gel(z,2)))
1701 113225 : return gerepileuptoint((pari_sp)(z+3), gel(z,1));
1702 : }
1703 : else
1704 : { /* naive 4M formula: avoid all loss of accuracy */
1705 254373285 : p1 = gmul(xr,yr);
1706 254334448 : p2 = gmul(xi,yi);
1707 254337682 : p3 = gmul(xr,yi);
1708 254341522 : p4 = gmul(xi,yr);
1709 254342053 : tetpil = avma;
1710 254342053 : gel(z,1) = gsub(p1,p2);
1711 254209266 : gel(z,2) = gadd(p3,p4);
1712 254216104 : if (isintzero(gel(z,2)))
1713 : {
1714 50575 : cgiv(gel(z,2));
1715 50575 : return gerepileupto((pari_sp)(z+3), gel(z,1));
1716 : }
1717 : }
1718 : #endif
1719 :
1720 254602755 : gerepilecoeffssp(av,tetpil, z+1,2); return z;
1721 : }
1722 : /* x,y PADIC */
1723 : static GEN
1724 17502482 : mulpp(GEN x, GEN y) {
1725 17502482 : long l = valp(x) + valp(y);
1726 : pari_sp av;
1727 : GEN z, t;
1728 17502482 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
1729 17502456 : if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
1730 17275788 : if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
1731 :
1732 16929230 : t = (precp(x) > precp(y))? y: x;
1733 16929230 : z = cgetp(t); setvalp(z,l); av = avma;
1734 16929461 : affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
1735 16929360 : set_avma(av); return z;
1736 : }
1737 : /* x,y QUAD */
1738 : static GEN
1739 1106 : mulqq(GEN x, GEN y) {
1740 1106 : GEN z = cgetg(4,t_QUAD);
1741 1106 : GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
1742 : pari_sp av, tetpil;
1743 1106 : if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
1744 :
1745 1106 : gel(z,1) = ZX_copy(P); av = avma;
1746 1106 : p2 = gmul(gel(x,2),gel(y,2));
1747 1106 : p3 = gmul(gel(x,3),gel(y,3));
1748 1106 : p1 = gmul(gneg_i(c),p3);
1749 :
1750 1106 : if (signe(b))
1751 987 : p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
1752 : else
1753 : {
1754 119 : p3 = gmul(gel(x,2),gel(y,3));
1755 119 : p4 = gmul(gel(x,3),gel(y,2));
1756 : }
1757 1106 : tetpil = avma;
1758 1106 : gel(z,2) = gadd(p2,p1);
1759 1106 : gel(z,3) = gadd(p4,p3);
1760 1106 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
1761 : }
1762 :
1763 : GEN
1764 15626814 : mulcxI(GEN x)
1765 : {
1766 : GEN z;
1767 15626814 : switch(typ(x))
1768 : {
1769 2627889 : case t_INT: case t_REAL: case t_FRAC:
1770 2627889 : return mkcomplex(gen_0, x);
1771 12999477 : case t_COMPLEX:
1772 12999477 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1773 12986036 : z = cgetg(3,t_COMPLEX);
1774 12986265 : gel(z,1) = gneg(gel(x,2));
1775 12986811 : gel(z,2) = gel(x,1); return z;
1776 48 : default:
1777 48 : return gmul(gen_I(), x);
1778 : }
1779 : }
1780 : GEN
1781 3079334 : mulcxmI(GEN x)
1782 : {
1783 : GEN z;
1784 3079334 : switch(typ(x))
1785 : {
1786 116080 : case t_INT: case t_REAL: case t_FRAC:
1787 116080 : return mkcomplex(gen_0, gneg(x));
1788 2760350 : case t_COMPLEX:
1789 2760350 : if (isintzero(gel(x,1))) return gel(x,2);
1790 2722165 : z = cgetg(3,t_COMPLEX);
1791 2722140 : gel(z,1) = gel(x,2);
1792 2722140 : gel(z,2) = gneg(gel(x,1)); return z;
1793 202904 : default:
1794 202904 : return gmul(mkcomplex(gen_0, gen_m1), x);
1795 : }
1796 : }
1797 : /* x * I^k */
1798 : GEN
1799 5523052 : mulcxpowIs(GEN x, long k)
1800 : {
1801 5523052 : switch (k & 3)
1802 : {
1803 1369089 : case 1: return mulcxI(x);
1804 1362094 : case 2: return gneg(x);
1805 1363687 : case 3: return mulcxmI(x);
1806 : }
1807 1428182 : return x;
1808 : }
1809 :
1810 : /* caller will assume l > 2 later */
1811 : static GEN
1812 5706235 : init_ser(long l, long v, long e)
1813 : {
1814 5706235 : GEN z = cgetg(l, t_SER);
1815 5706235 : z[1] = evalvalser(e) | evalvarn(v) | evalsigne(1); return z;
1816 : }
1817 :
1818 : /* fill in coefficients of t_SER z from coeffs of t_POL y */
1819 : static GEN
1820 5693888 : fill_ser(GEN z, GEN y)
1821 : {
1822 5693888 : long i, lx = lg(z), ly = lg(y); /* lx > 2 */
1823 5693888 : if (ly >= lx) {
1824 25570501 : for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
1825 : } else {
1826 335415 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1827 239573 : for ( ; i < lx; i++) gel(z,i) = gen_0;
1828 : }
1829 : /* dangerous case (already normalized), don't use normalizeser */
1830 5693888 : if (ly == 3 && !signe(y)) { setsigne(z, 0); return z; }
1831 5693790 : return normalizeser(z);
1832 : }
1833 : /* assume typ(x) = t_VEC */
1834 : static int
1835 98 : is_ext_qfr(GEN x)
1836 91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
1837 189 : && typ(gel(x,2)) == t_REAL; }
1838 :
1839 : GEN
1840 8517887033 : gmul(GEN x, GEN y)
1841 : {
1842 : long tx, ty, lx, ly, vx, vy, i, l;
1843 : pari_sp av, tetpil;
1844 : GEN z, p1;
1845 :
1846 8517887033 : if (x == y) return gsqr(x);
1847 7639709350 : tx = typ(x); ty = typ(y);
1848 7639709350 : if (tx == ty) switch(tx)
1849 : {
1850 3623166036 : case t_INT: return mulii(x,y);
1851 1867800129 : case t_REAL: return mulrr(x,y);
1852 1780475 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1853 1780475 : z = cgetg(3,t_INTMOD);
1854 1780475 : if (X==Y || equalii(X,Y))
1855 1780440 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1856 35 : gel(z,1) = gcdii(X,Y);
1857 35 : warn_coercion(X,Y,gel(z,1));
1858 35 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1859 35 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
1860 : }
1861 23418660 : case t_FRAC:
1862 : {
1863 23418660 : GEN x1 = gel(x,1), x2 = gel(x,2);
1864 23418660 : GEN y1 = gel(y,1), y2 = gel(y,2);
1865 23418660 : z=cgetg(3,t_FRAC);
1866 23418662 : p1 = gcdii(x1, y2);
1867 23418662 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1868 23418662 : p1 = gcdii(x2, y1);
1869 23418661 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1870 23418660 : tetpil = avma;
1871 23418660 : gel(z,2) = mulii(x2,y2);
1872 23418659 : gel(z,1) = mulii(x1,y1);
1873 23418661 : fix_frac_if_int_GC(z,tetpil); return z;
1874 : }
1875 264753581 : case t_COMPLEX: return mulcc(x, y);
1876 12510935 : case t_PADIC: return mulpp(x, y);
1877 882 : case t_QUAD: return mulqq(x, y);
1878 14545974 : case t_FFELT: return FF_mul(x, y);
1879 11003949 : case t_POLMOD:
1880 11003949 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1881 8328843 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1882 2675106 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1883 48071411 : case t_POL:
1884 48071411 : vx = varn(x);
1885 48071411 : vy = varn(y);
1886 48071411 : if (vx != vy) {
1887 4854818 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1888 2071778 : else return RgX_Rg_mul(y, x);
1889 : }
1890 43216593 : return RgX_mul(x, y);
1891 :
1892 4411812 : case t_SER: {
1893 4411812 : vx = varn(x);
1894 4411812 : vy = varn(y);
1895 4411812 : if (vx != vy) {
1896 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1897 3675 : return mul_ser_scal(y, x);
1898 : }
1899 4408137 : lx = minss(lg(x), lg(y));
1900 4408137 : if (lx == 2) return zeroser(vx, valser(x)+valser(y));
1901 4408123 : av = avma; z = init_ser(lx, vx, valser(x)+valser(y));
1902 4408123 : x = ser2pol_i(x, lx);
1903 4408123 : y = ser2pol_i(y, lx);
1904 4408123 : y = RgXn_mul(x, y, lx-2);
1905 4408123 : return gerepilecopy(av, fill_ser(z,y));
1906 : }
1907 42 : case t_VEC:
1908 42 : if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
1909 : /* fall through, handle extended t_QFB */
1910 47804 : case t_QFB: return qfbcomp(x,y);
1911 6720535 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1912 4114630 : case t_MAT: return RgM_mul(x, y);
1913 :
1914 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1915 1631 : z = cgetg_copy(x, &l);
1916 1631 : if (l != lg(y)) break;
1917 17325 : for (i=1; i<l; i++)
1918 : {
1919 15694 : long yi = y[i];
1920 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1921 15694 : z[i] = x[yi];
1922 : }
1923 1631 : return z;
1924 :
1925 0 : default:
1926 0 : pari_err_TYPE2("*",x,y);
1927 : }
1928 : /* tx != ty */
1929 1759875946 : if (is_const_t(ty) && is_const_t(tx)) {
1930 1626479125 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1931 1626479125 : switch(tx) {
1932 1497658323 : case t_INT:
1933 : switch(ty)
1934 : {
1935 987260740 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1936 1224341 : case t_INTMOD:
1937 1224341 : z = cgetg(3, t_INTMOD);
1938 1224339 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1939 72545009 : case t_FRAC:
1940 72545009 : if (!signe(x)) return gen_0;
1941 55243606 : z=cgetg(3,t_FRAC);
1942 55243803 : p1 = gcdii(x,gel(y,2));
1943 55243003 : if (equali1(p1))
1944 : {
1945 28963448 : set_avma((pari_sp)z);
1946 28963439 : gel(z,2) = icopy(gel(y,2));
1947 28963474 : gel(z,1) = mulii(gel(y,1), x);
1948 : }
1949 : else
1950 : {
1951 26279800 : x = diviiexact(x,p1); tetpil = avma;
1952 26279203 : gel(z,2) = diviiexact(gel(y,2), p1);
1953 26278994 : gel(z,1) = mulii(gel(y,1), x);
1954 26279507 : fix_frac_if_int_GC(z,tetpil);
1955 : }
1956 55243520 : return z;
1957 433892547 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1958 4196328 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1959 1701 : case t_QUAD: return mulRq(x,y);
1960 1627576 : case t_FFELT: return FF_Z_mul(y,x);
1961 : }
1962 :
1963 : case t_REAL:
1964 121130402 : switch(ty)
1965 : {
1966 37507656 : case t_FRAC: return mulrfrac(x, y);
1967 83622704 : case t_COMPLEX: return mulRc(x, y);
1968 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1969 21 : default: pari_err_TYPE2("*",x,y);
1970 : }
1971 :
1972 8022 : case t_INTMOD:
1973 : switch(ty)
1974 : {
1975 7133 : case t_FRAC: { GEN X = gel(x,1);
1976 7133 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1977 7133 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1978 : }
1979 49 : case t_COMPLEX: return mulRc_direct(x,y);
1980 427 : case t_PADIC: { GEN X = gel(x,1);
1981 427 : z = cgetg(3, t_INTMOD);
1982 427 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1983 : }
1984 63 : case t_QUAD: return mulRq(x, y);
1985 350 : case t_FFELT:
1986 350 : if (!equalii(gel(x,1),FF_p_i(y)))
1987 0 : pari_err_OP("*",x,y);
1988 350 : return FF_Z_mul(y,gel(x,2));
1989 : }
1990 :
1991 : case t_FRAC:
1992 : switch(ty)
1993 : {
1994 5235229 : case t_COMPLEX: return mulRc(x, y);
1995 2582763 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
1996 343 : case t_QUAD: return mulRq(x, y);
1997 2051 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
1998 : }
1999 :
2000 : case t_FFELT:
2001 35 : pari_err_TYPE2("*",x,y);
2002 :
2003 21 : case t_COMPLEX:
2004 : switch(ty)
2005 : {
2006 14 : case t_PADIC:
2007 14 : return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
2008 7 : case t_QUAD:
2009 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
2010 7 : return mulqf(y, x, lx);
2011 : }
2012 :
2013 : case t_PADIC: /* ty == t_QUAD */
2014 28 : return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
2015 : }
2016 : }
2017 :
2018 135604576 : if (is_matvec_t(ty)) switch(tx)
2019 : {
2020 122599 : case t_VEC:
2021 : switch(ty) {
2022 15540 : case t_COL: return RgV_RgC_mul(x,y);
2023 107059 : case t_MAT: return RgV_RgM_mul(x,y);
2024 : }
2025 0 : break;
2026 1687 : case t_COL:
2027 : switch(ty) {
2028 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2029 0 : case t_MAT: return RgC_RgM_mul(x,y);
2030 : }
2031 0 : break;
2032 825137 : case t_MAT:
2033 : switch(ty) {
2034 0 : case t_VEC: return RgM_RgV_mul(x,y);
2035 825137 : case t_COL: return RgM_RgC_mul(x,y);
2036 : }
2037 : default:
2038 7457802 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
2039 7457826 : return mul_self_scal(y, x);
2040 : }
2041 130588353 : if (is_matvec_t(tx))
2042 : {
2043 3084297 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2044 3084715 : return mul_self_scal(x, y);
2045 : }
2046 127504627 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2047 : /* tx < ty, !ismatvec(x and y) */
2048 :
2049 127504627 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2050 3403942 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2051 124100685 : if (is_scalar_t(tx)) {
2052 120809312 : if (tx == t_POLMOD) {
2053 3106910 : vx = varn(gel(x,1));
2054 3106910 : vy = gvar(y);
2055 3106910 : if (vx != vy) {
2056 2872620 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2057 23324 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2058 : }
2059 : /* error if ty == t_SER */
2060 234290 : av = avma; y = gmod(y, gel(x,1));
2061 234283 : return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2062 : }
2063 117702402 : return mul_scal(y, x, ty);
2064 : }
2065 :
2066 : /* x and y are not scalars, nor matvec */
2067 3291256 : vx = gvar(x);
2068 3291286 : vy = gvar(y);
2069 3291286 : if (vx != vy) /* x or y is treated as a scalar */
2070 2780439 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2071 2780439 : : mul_scal(y, x, ty);
2072 : /* vx = vy */
2073 1872272 : switch(tx)
2074 : {
2075 1872237 : case t_POL:
2076 : switch (ty)
2077 : {
2078 6741 : case t_SER:
2079 : {
2080 : long v;
2081 6741 : av = avma; v = RgX_valrem(x, &x);
2082 6741 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2083 6727 : v += valser(y); ly = lg(y);
2084 6727 : if (ly == 2) { set_avma(av); return zeroser(vx, v); }
2085 6496 : if (degpol(x))
2086 : {
2087 2030 : z = init_ser(ly, vx, v);
2088 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2089 2030 : return gerepilecopy(av, fill_ser(z, x));
2090 : }
2091 : /* take advantage of x = c*t^v */
2092 4466 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2093 4466 : setvalser(y, v); return y;
2094 : }
2095 :
2096 1865496 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2097 : }
2098 0 : break;
2099 :
2100 35 : case t_SER:
2101 : switch (ty)
2102 : {
2103 35 : case t_RFRAC:
2104 35 : av = avma;
2105 35 : return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2106 : }
2107 0 : break;
2108 : }
2109 0 : pari_err_TYPE2("*",x,y);
2110 : return NULL; /* LCOV_EXCL_LINE */
2111 : }
2112 :
2113 : /* return a nonnormalized result */
2114 : GEN
2115 112731 : sqr_ser_part(GEN x, long l1, long l2)
2116 : {
2117 : long i, j, l;
2118 : pari_sp av;
2119 : GEN Z, z, p1, p2;
2120 : long mi;
2121 112731 : if (l2 < l1) return zeroser(varn(x), 2*valser(x));
2122 112717 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2123 112717 : Z = cgetg(l2-l1+3, t_SER);
2124 112717 : Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
2125 112717 : z = Z + 2-l1;
2126 112717 : x += 2; mi = 0;
2127 425390 : for (i=0; i<l1; i++)
2128 : {
2129 312673 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2130 : }
2131 :
2132 720212 : for (i=l1; i<=l2; i++)
2133 : {
2134 607495 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2135 607495 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2136 1256134 : for (j=i-mi; j<=minss(l,mi); j++)
2137 648639 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2138 607495 : p1 = gshift(p1,1);
2139 607495 : if ((i&1) == 0 && p2[i>>1])
2140 93585 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2141 607495 : gel(z,i) = gerepileupto(av,p1);
2142 : }
2143 112717 : return Z;
2144 : }
2145 :
2146 : GEN
2147 1122280252 : gsqr(GEN x)
2148 : {
2149 : long i, lx;
2150 : pari_sp av, tetpil;
2151 : GEN z, p1, p2, p3, p4;
2152 :
2153 1122280252 : switch(typ(x))
2154 : {
2155 909768827 : case t_INT: return sqri(x);
2156 195242574 : case t_REAL: return sqrr(x);
2157 142490 : case t_INTMOD: { GEN X = gel(x,1);
2158 142490 : z = cgetg(3,t_INTMOD);
2159 142491 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
2160 142490 : gel(z,1) = icopy(X); return z;
2161 : }
2162 4141136 : case t_FRAC: return sqrfrac(x);
2163 :
2164 8109074 : case t_COMPLEX:
2165 8109074 : if (isintzero(gel(x,1))) {
2166 238822 : av = avma;
2167 238822 : return gerepileupto(av, gneg(gsqr(gel(x,2))));
2168 : }
2169 7870259 : z = cgetg(3,t_COMPLEX); av = avma;
2170 7870239 : p1 = gadd(gel(x,1),gel(x,2));
2171 7870052 : p2 = gsub(gel(x,1), gel(x,2));
2172 7870059 : p3 = gmul(gel(x,1),gel(x,2));
2173 7870200 : tetpil = avma;
2174 7870200 : gel(z,1) = gmul(p1,p2);
2175 7870213 : gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
2176 :
2177 4746 : case t_PADIC:
2178 4746 : z = cgetg(5,t_PADIC);
2179 4746 : i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
2180 4746 : if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
2181 4746 : z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
2182 4746 : gel(z,2) = icopy(gel(x,2));
2183 4746 : gel(z,3) = shifti(gel(x,3), i); av = avma;
2184 4746 : gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
2185 4746 : return z;
2186 :
2187 70 : case t_QUAD: z = cgetg(4,t_QUAD);
2188 70 : p1 = gel(x,1);
2189 70 : gel(z,1) = ZX_copy(p1); av = avma;
2190 70 : p2 = gsqr(gel(x,2));
2191 70 : p3 = gsqr(gel(x,3));
2192 70 : p4 = gmul(gneg_i(gel(p1,2)),p3);
2193 :
2194 70 : if (gequal0(gel(p1,3)))
2195 : {
2196 7 : tetpil = avma;
2197 7 : gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
2198 7 : av = avma;
2199 7 : p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
2200 7 : gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
2201 : }
2202 :
2203 63 : p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
2204 63 : tetpil = avma;
2205 63 : gel(z,2) = gadd(p2,p4);
2206 63 : gel(z,3) = gadd(p1,p3);
2207 63 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
2208 :
2209 484415 : case t_POLMOD:
2210 484415 : return sqr_polmod(gel(x,1), gel(x,2));
2211 :
2212 2956601 : case t_FFELT: return FF_sqr(x);
2213 :
2214 1327856 : case t_POL: return RgX_sqr(x);
2215 :
2216 35014 : case t_SER:
2217 35014 : lx = lg(x);
2218 35014 : if (ser_isexactzero(x)) {
2219 21 : GEN z = gcopy(x);
2220 21 : setvalser(z, 2*valser(x));
2221 21 : return z;
2222 : }
2223 34993 : if (lx < 40)
2224 34734 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2225 : else
2226 : {
2227 259 : pari_sp av = avma;
2228 259 : GEN z = init_ser(lx, varn(x), 2*valser(x));
2229 259 : x = ser2pol_i(x, lx);
2230 259 : x = RgXn_sqr(x, lx-2);
2231 259 : return gerepilecopy(av, fill_ser(z,x));
2232 : }
2233 :
2234 70 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2235 70 : gel(z,1) = gsqr(gel(x,1));
2236 70 : gel(z,2) = gsqr(gel(x,2)); return z;
2237 :
2238 1113 : case t_MAT: return RgM_sqr(x);
2239 28 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
2240 : /* fall through handle extended t_QFB */
2241 81919 : case t_QFB: return qfbsqr(x);
2242 658 : case t_VECSMALL:
2243 658 : z = cgetg_copy(x, &lx);
2244 16289 : for (i=1; i<lx; i++)
2245 : {
2246 15631 : long xi = x[i];
2247 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2248 15631 : z[i] = x[xi];
2249 : }
2250 658 : return z;
2251 : }
2252 7 : pari_err_TYPE2("*",x,x);
2253 : return NULL; /* LCOV_EXCL_LINE */
2254 : }
2255 :
2256 : /********************************************************************/
2257 : /** **/
2258 : /** DIVISION **/
2259 : /** **/
2260 : /********************************************************************/
2261 : static GEN
2262 13169 : div_rfrac_scal(GEN x, GEN y)
2263 : {
2264 13169 : pari_sp av = avma;
2265 13169 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2266 13169 : return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
2267 : }
2268 : static GEN
2269 37515 : div_scal_rfrac(GEN x, GEN y)
2270 : {
2271 37515 : GEN y1 = gel(y,1), y2 = gel(y,2);
2272 37515 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2273 : {
2274 14 : if (degpol(y1))
2275 : {
2276 14 : pari_sp av = avma;
2277 14 : GEN _1 = Rg_get_1(x);
2278 14 : if (transtype(_1)) y1 = gmul(y1, _1);
2279 14 : return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
2280 : }
2281 0 : y1 = gel(y1,2);
2282 : }
2283 37501 : return RgX_Rg_mul(y2, gdiv(x,y1));
2284 : }
2285 : static GEN
2286 1186775 : div_rfrac(GEN x, GEN y)
2287 1186775 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2288 :
2289 : /* x != 0 */
2290 : static GEN
2291 1338593 : div_ser_scal(GEN x, GEN t)
2292 : {
2293 1338593 : if (ser_isexactzero(x))
2294 : {
2295 28 : GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gdiv(gel(x,2), t), varn(x), 1);
2296 28 : setvalser(z, valser(x)); return z;
2297 : }
2298 6202113 : pari_APPLY_ser(gdiv(gel(x,i), t));
2299 : }
2300 : GEN
2301 7630 : ser_normalize(GEN x)
2302 : {
2303 7630 : long i, lx = lg(x);
2304 : GEN c, z;
2305 7630 : if (lx == 2) return x;
2306 7630 : c = gel(x,2); if (gequal1(c)) return x;
2307 7553 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2308 108528 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2309 7553 : return z;
2310 : }
2311 :
2312 : /* y != 0 */
2313 : static GEN
2314 6685203 : div_T_scal(GEN x, GEN y, long tx) {
2315 6685203 : switch(tx)
2316 : {
2317 5335520 : case t_POL: return RgX_Rg_div(x, y);
2318 1338586 : case t_SER: return div_ser_scal(x, y);
2319 11097 : case t_RFRAC: return div_rfrac_scal(x,y);
2320 : }
2321 0 : pari_err_TYPE2("/",x,y);
2322 : return NULL; /* LCOV_EXCL_LINE */
2323 : }
2324 :
2325 : static GEN
2326 9257890 : div_scal_pol(GEN x, GEN y) {
2327 9257890 : long ly = lg(y);
2328 : pari_sp av;
2329 : GEN _1;
2330 9257890 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2331 9180196 : if (isrationalzero(x)) return zeropol(varn(y));
2332 7130228 : av = avma;
2333 7130228 : _1 = Rg_get_1(x); if (transtype(_1)) y = gmul(y, _1);
2334 7130228 : return gerepileupto(av, gred_rfrac_simple(x,y));
2335 : }
2336 : static GEN
2337 18550 : div_scal_ser(GEN x, GEN y)
2338 : {
2339 18550 : pari_sp av = avma;
2340 18550 : GEN _1 = Rg_get_1(x);
2341 18550 : if (transtype(_1)) y = gmul(y, _1);
2342 18550 : return gerepileupto(av, gmul(x, ser_inv(y)));
2343 : }
2344 : static GEN
2345 9267125 : div_scal_T(GEN x, GEN y, long ty) {
2346 9267125 : switch(ty)
2347 : {
2348 9211823 : case t_POL: return div_scal_pol(x, y);
2349 18543 : case t_SER: return div_scal_ser(x, y);
2350 36759 : case t_RFRAC: return div_scal_rfrac(x, y);
2351 : }
2352 0 : pari_err_TYPE2("/",x,y);
2353 : return NULL; /* LCOV_EXCL_LINE */
2354 : }
2355 :
2356 : /* assume tx = ty = t_SER, same variable vx */
2357 : static GEN
2358 771510 : div_ser(GEN x, GEN y, long vx)
2359 : {
2360 771510 : long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
2361 771510 : GEN y0 = y, z;
2362 771510 : pari_sp av = avma;
2363 :
2364 771510 : if (!signe(y)) pari_err_INV("div_ser", y);
2365 771503 : if (ser_isexactzero(x))
2366 : {
2367 59892 : if (lx == 2) return zeroser(vx, v);
2368 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2369 147 : setvalser(z, v); return z;
2370 : }
2371 711611 : if (lx < ly) ly = lx;
2372 711611 : y = ser2pol_i_normalize(y, ly, &e);
2373 711611 : if (e)
2374 : {
2375 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2376 7 : ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
2377 : }
2378 711611 : z = init_ser(ly, vx, v);
2379 711611 : if (ly == 3)
2380 : {
2381 12347 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2382 12347 : if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
2383 12347 : return gerepileupto(av, z);
2384 : }
2385 699264 : x = ser2pol_i(x, ly);
2386 699264 : y = RgXn_div_i(x, y, ly-2);
2387 699264 : return gerepilecopy(av, fill_ser(z,y));
2388 : }
2389 : /* x,y compatible PADIC */
2390 : static GEN
2391 13212908 : divpp(GEN x, GEN y) {
2392 : pari_sp av;
2393 : long a, b;
2394 : GEN z, M;
2395 :
2396 13212908 : if (!signe(gel(y,4))) pari_err_INV("divpp",y);
2397 13212907 : if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
2398 13212578 : a = precp(x);
2399 13212578 : b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
2400 13212578 : z = cgetg(5, t_PADIC);
2401 13212242 : z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
2402 13212167 : gel(z,2) = icopy(gel(x,2));
2403 13211853 : gel(z,3) = icopy(M); av = avma;
2404 13211442 : gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
2405 13212046 : return z;
2406 : }
2407 : static GEN
2408 36764 : div_polmod_same(GEN T, GEN x, GEN y)
2409 : {
2410 36764 : long v = varn(T);
2411 36764 : GEN a, z = cgetg(3, t_POLMOD);
2412 36764 : gel(z,1) = RgX_copy(T);
2413 36764 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2414 11130 : a = gdiv(x, y);
2415 25634 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2416 119 : {
2417 119 : pari_sp av = avma;
2418 119 : a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
2419 : }
2420 25515 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2421 5271 : {
2422 5271 : pari_sp av = avma;
2423 5271 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2424 5271 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2425 5271 : a = gerepileupto(av, a);
2426 : }
2427 : else
2428 : {
2429 20244 : pari_sp av = avma;
2430 20244 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2431 20244 : a = gerepileupto(av, a);
2432 : }
2433 36764 : gel(z,2) = a; return z;
2434 : }
2435 : GEN
2436 432525457 : gdiv(GEN x, GEN y)
2437 : {
2438 432525457 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2439 : pari_sp av, tetpil;
2440 : GEN z, p1, p2;
2441 :
2442 432525457 : if (tx == ty) switch(tx)
2443 : {
2444 88660239 : case t_INT:
2445 88660239 : return Qdivii(x,y);
2446 :
2447 93620588 : case t_REAL: return divrr(x,y);
2448 25746 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2449 25746 : z = cgetg(3,t_INTMOD);
2450 25746 : if (X==Y || equalii(X,Y))
2451 25732 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2452 14 : gel(z,1) = gcdii(X,Y);
2453 14 : warn_coercion(X,Y,gel(z,1));
2454 14 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2455 14 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
2456 : }
2457 4041449 : case t_FRAC: {
2458 4041449 : GEN x1 = gel(x,1), x2 = gel(x,2);
2459 4041449 : GEN y1 = gel(y,1), y2 = gel(y,2);
2460 4041449 : z = cgetg(3, t_FRAC);
2461 4041449 : p1 = gcdii(x1, y1);
2462 4041448 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2463 4041448 : p1 = gcdii(x2, y2);
2464 4041447 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2465 4041447 : tetpil = avma;
2466 4041447 : gel(z,2) = mulii(x2,y1);
2467 4041446 : gel(z,1) = mulii(x1,y2);
2468 4041445 : normalize_frac(z);
2469 4041445 : fix_frac_if_int_GC(z,tetpil);
2470 4041450 : return z;
2471 : }
2472 9200956 : case t_COMPLEX:
2473 9200956 : if (isintzero(gel(y,1)))
2474 : {
2475 198044 : y = gel(y,2);
2476 198044 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2477 121135 : z = cgetg(3,t_COMPLEX);
2478 121135 : gel(z,1) = gdiv(gel(x,2), y);
2479 121135 : av = avma;
2480 121135 : gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
2481 121135 : return z;
2482 : }
2483 9002907 : av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
2484 9002894 : return gerepile(av, tetpil, gdiv(p2,p1));
2485 :
2486 587420 : case t_PADIC:
2487 587420 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
2488 587420 : return divpp(x, y);
2489 :
2490 238 : case t_QUAD:
2491 238 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2492 224 : av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
2493 224 : return gerepile(av, tetpil, gdiv(p2,p1));
2494 :
2495 243520 : case t_FFELT: return FF_div(x,y);
2496 :
2497 41139 : case t_POLMOD:
2498 41139 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2499 36764 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2500 : else {
2501 4375 : av = avma;
2502 4375 : z = gerepileupto(av, gmul(x, ginv(y)));
2503 : }
2504 41139 : return z;
2505 :
2506 18453300 : case t_POL:
2507 18453300 : vx = varn(x);
2508 18453300 : vy = varn(y);
2509 18453300 : if (vx != vy) {
2510 102172 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2511 46067 : else return div_scal_pol(x, y);
2512 : }
2513 18351128 : if (!signe(y)) pari_err_INV("gdiv",y);
2514 18351128 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2515 18228836 : av = avma;
2516 18228836 : return gerepileupto(av, gred_rfrac2(x,y));
2517 :
2518 771531 : case t_SER:
2519 771531 : vx = varn(x);
2520 771531 : vy = varn(y);
2521 771531 : if (vx != vy) {
2522 21 : if (varncmp(vx, vy) < 0)
2523 : {
2524 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2525 7 : return div_ser_scal(x, y);
2526 : }
2527 7 : return div_scal_ser(x, y);
2528 : }
2529 771510 : return div_ser(x, y, vx);
2530 1189463 : case t_RFRAC:
2531 1189463 : vx = varn(gel(x,2));
2532 1189463 : vy = varn(gel(y,2));
2533 1189463 : if (vx != vy) {
2534 2688 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2535 756 : else return div_scal_rfrac(x, y);
2536 : }
2537 1186775 : return div_rfrac(x,y);
2538 :
2539 21 : case t_VEC: /* handle extended t_QFB */
2540 21 : case t_QFB: av = avma; return gerepileupto(av, qfbcomp(x, ginv(y)));
2541 :
2542 28 : case t_MAT:
2543 28 : p1 = RgM_div(x,y);
2544 28 : if (!p1) pari_err_INV("gdiv",y);
2545 21 : return p1;
2546 :
2547 0 : default: pari_err_TYPE2("/",x,y);
2548 : }
2549 :
2550 215694664 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2551 : {
2552 3777812 : long s = signe(x);
2553 3777812 : if (!s) {
2554 520742 : if (gequal0(y)) pari_err_INV("gdiv",y);
2555 520742 : switch (ty)
2556 : {
2557 517578 : default: return gen_0;
2558 28 : case t_INTMOD:
2559 28 : z = cgetg(3,t_INTMOD);
2560 28 : gel(z,1) = icopy(gel(y,1));
2561 28 : gel(z,2) = gen_0; return z;
2562 3136 : case t_FFELT: return FF_zero(y);
2563 : }
2564 : }
2565 3257070 : if (is_pm1(x)) {
2566 1307700 : if (s > 0) return ginv(y);
2567 221065 : av = avma; return gerepileupto(av, ginv(gneg(y)));
2568 : }
2569 1949370 : switch(ty)
2570 : {
2571 639063 : case t_REAL: return divir(x,y);
2572 21 : case t_INTMOD:
2573 21 : z = cgetg(3, t_INTMOD);
2574 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2575 794565 : case t_FRAC:
2576 794565 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2577 794565 : if (equali1(p1))
2578 : {
2579 273786 : set_avma((pari_sp)z);
2580 273786 : gel(z,2) = icopy(gel(y,1));
2581 273786 : gel(z,1) = mulii(gel(y,2), x);
2582 273786 : normalize_frac(z);
2583 273786 : fix_frac_if_int(z);
2584 : }
2585 : else
2586 : {
2587 520779 : x = diviiexact(x,p1); tetpil = avma;
2588 520779 : gel(z,2) = diviiexact(gel(y,1), p1);
2589 520779 : gel(z,1) = mulii(gel(y,2), x);
2590 520779 : normalize_frac(z);
2591 520779 : fix_frac_if_int_GC(z,tetpil);
2592 : }
2593 794565 : return z;
2594 :
2595 462 : case t_FFELT: return Z_FF_div(x,y);
2596 513526 : case t_COMPLEX: return divRc(x,y);
2597 1505 : case t_PADIC: return divTp(x, y);
2598 231 : case t_QUAD:
2599 231 : av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
2600 231 : return gerepile(av, tetpil, gdiv(p2,p1));
2601 : }
2602 : }
2603 211916849 : if (gequal0(y))
2604 : {
2605 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2606 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2607 : }
2608 :
2609 211926958 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2610 : {
2611 172564085 : case t_REAL:
2612 172564085 : switch(ty)
2613 : {
2614 170203319 : case t_INT: return divri(x,y);
2615 522503 : case t_FRAC:
2616 522503 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2617 522503 : return gerepileuptoleaf(av, z);
2618 1838293 : case t_COMPLEX: return divRc(x, y);
2619 42 : case t_QUAD: return divfq(x, y, realprec(x));
2620 18 : default: pari_err_TYPE2("/",x,y);
2621 : }
2622 :
2623 280 : case t_INTMOD:
2624 : switch(ty)
2625 : {
2626 196 : case t_INT:
2627 196 : z = cgetg(3, t_INTMOD);
2628 196 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2629 28 : case t_FRAC: { GEN X = gel(x,1);
2630 28 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2631 28 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2632 : }
2633 7 : case t_FFELT:
2634 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2635 0 : pari_err_OP("/",x,y);
2636 7 : return Z_FF_div(gel(x,2),y);
2637 :
2638 0 : case t_COMPLEX:
2639 0 : av = avma;
2640 0 : return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
2641 :
2642 14 : case t_QUAD:
2643 14 : av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
2644 14 : return gerepile(av,tetpil, gdiv(p2,p1));
2645 :
2646 7 : case t_PADIC: { GEN X = gel(x,1);
2647 7 : z = cgetg(3, t_INTMOD);
2648 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2649 : }
2650 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2651 : }
2652 :
2653 : case t_FRAC:
2654 : switch(ty)
2655 : {
2656 2097616 : case t_INT: z = cgetg(3, t_FRAC);
2657 2097616 : p1 = gcdii(y,gel(x,1));
2658 2097616 : if (equali1(p1))
2659 : {
2660 828157 : set_avma((pari_sp)z); tetpil = 0;
2661 828157 : gel(z,1) = icopy(gel(x,1));
2662 : }
2663 : else
2664 : {
2665 1269459 : y = diviiexact(y,p1); tetpil = avma;
2666 1269459 : gel(z,1) = diviiexact(gel(x,1), p1);
2667 : }
2668 2097616 : gel(z,2) = mulii(gel(x,2),y);
2669 2097616 : normalize_frac(z);
2670 2097616 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2671 2097616 : return z;
2672 :
2673 59516 : case t_REAL:
2674 59516 : av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
2675 59516 : return gerepile(av, tetpil, divir(gel(x,1), p1));
2676 :
2677 7 : case t_INTMOD: { GEN Y = gel(y,1);
2678 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2679 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2680 : }
2681 :
2682 28 : case t_FFELT: av=avma;
2683 28 : return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2684 :
2685 861 : case t_COMPLEX: return divRc(x, y);
2686 :
2687 2141 : case t_PADIC:
2688 2141 : if (!signe(gel(x,1))) return gen_0;
2689 2141 : return divTp(x, y);
2690 :
2691 42 : case t_QUAD:
2692 42 : av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
2693 42 : return gerepile(av,tetpil,gdiv(p2,p1));
2694 : }
2695 :
2696 : case t_FFELT:
2697 161 : switch (ty)
2698 : {
2699 77 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2700 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2701 7 : case t_INTMOD:
2702 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2703 0 : pari_err_OP("/",x,y);
2704 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2705 49 : default:
2706 49 : pari_err_TYPE2("/",x,y);
2707 : }
2708 0 : break;
2709 :
2710 13507830 : case t_COMPLEX:
2711 : switch(ty)
2712 : {
2713 13507831 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2714 0 : case t_INTMOD: return mulRc_direct(ginv(y), x);
2715 0 : case t_PADIC:
2716 0 : return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
2717 0 : case t_QUAD:
2718 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2719 0 : return divfq(x, y, lx);
2720 : }
2721 :
2722 : case t_PADIC:
2723 : switch(ty)
2724 : {
2725 1205001 : case t_INT: case t_FRAC: { GEN p = gel(x,2);
2726 1204910 : return signe(gel(x,4))? divpT(x, y)
2727 2409911 : : zeropadic(p, valp(x) - Q_pval(y,p));
2728 : }
2729 7 : case t_INTMOD: { GEN Y = gel(y,1);
2730 7 : z = cgetg(3, t_INTMOD);
2731 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2732 : }
2733 14 : case t_COMPLEX: case t_QUAD:
2734 14 : av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
2735 14 : return gerepile(av,tetpil,gdiv(p1,p2));
2736 :
2737 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2738 : }
2739 :
2740 : case t_QUAD:
2741 : switch (ty)
2742 : {
2743 1190 : case t_INT: case t_INTMOD: case t_FRAC:
2744 1190 : z = cgetg(4,t_QUAD);
2745 1190 : gel(z,1) = ZX_copy(gel(x,1));
2746 1190 : gel(z,2) = gdiv(gel(x,2), y);
2747 1190 : gel(z,3) = gdiv(gel(x,3), y); return z;
2748 63 : case t_REAL: return divqf(x, y, realprec(y));
2749 14 : case t_PADIC: return divTp(x, y);
2750 0 : case t_COMPLEX:
2751 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2752 0 : return divqf(x, y, ly);
2753 : }
2754 : }
2755 22487899 : switch(ty) {
2756 581820 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2757 581820 : return gmul(x, ginv(y)); /* missing gerepile, for speed */
2758 42 : case t_MAT:
2759 42 : av = avma; p1 = RgM_inv(y);
2760 35 : if (!p1) pari_err_INV("gdiv",y);
2761 28 : return gerepileupto(av, gmul(x, p1));
2762 0 : case t_VEC: case t_COL:
2763 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2764 0 : pari_err_TYPE2("/",x,y);
2765 : }
2766 21906037 : switch(tx) {
2767 3843790 : case t_VEC: case t_COL: case t_MAT:
2768 12997121 : pari_APPLY_same(gdiv(gel(x,i),y));
2769 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2770 0 : pari_err_TYPE2("/",x,y);
2771 : }
2772 :
2773 18062247 : vy = gvar(y);
2774 18063251 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2775 209674 : vx = varn(X);
2776 209674 : if (vx != vy) {
2777 209100 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2778 208771 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2779 : }
2780 : /* y is POL, SER or RFRAC */
2781 574 : av = avma;
2782 574 : switch(ty)
2783 : {
2784 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2785 574 : default: y = ginvmod(gmod(y,X), X);
2786 : }
2787 567 : return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
2788 : }
2789 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2790 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2791 : * variable NO_VARIABLE, then the operation is incorrect */
2792 17853577 : vx = gvar(x);
2793 17853567 : if (vx != vy) { /* includes cases where one is scalar */
2794 15952006 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2795 9266795 : else return div_scal_T(x, y, ty);
2796 : }
2797 1901561 : switch(tx)
2798 : {
2799 1046213 : case t_POL:
2800 : switch(ty)
2801 : {
2802 28 : case t_SER:
2803 : {
2804 28 : GEN y0 = y;
2805 : long v;
2806 28 : av = avma; v = RgX_valrem(x, &x);
2807 28 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2808 14 : v -= valser(y); ly = lg(y); /* > 2 */
2809 14 : y = ser2pol_i_normalize(y, ly, &i);
2810 14 : if (i)
2811 : {
2812 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2813 7 : ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
2814 : }
2815 7 : z = init_ser(ly, vx, v);
2816 7 : return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, ly-2)));
2817 : }
2818 :
2819 1046185 : case t_RFRAC:
2820 : {
2821 1046185 : GEN y1 = gel(y,1), y2 = gel(y,2);
2822 1046185 : if (typ(y1) == t_POL && varn(y1) == vx)
2823 1171 : return mul_rfrac_scal(y2, y1, x);
2824 1045014 : av = avma;
2825 1045014 : return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2826 : }
2827 : }
2828 0 : break;
2829 :
2830 584324 : case t_SER:
2831 : switch(ty)
2832 : {
2833 584310 : case t_POL:
2834 : {
2835 584310 : long v = valser(x);
2836 584310 : lx = lg(x);
2837 584310 : if (lx == 2) return zeroser(vx, v - RgX_val(y));
2838 584205 : av = avma;
2839 584205 : x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
2840 584205 : z = init_ser(lx, vx, v);
2841 584205 : if (!signe(x)) setsigne(z,0);
2842 584205 : return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
2843 : }
2844 14 : case t_RFRAC:
2845 14 : av = avma;
2846 14 : return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2847 : }
2848 0 : break;
2849 :
2850 271003 : case t_RFRAC:
2851 : switch(ty)
2852 : {
2853 271003 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2854 0 : case t_SER:
2855 0 : av = avma;
2856 0 : return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2857 : }
2858 0 : break;
2859 : }
2860 21 : pari_err_TYPE2("/",x,y);
2861 : return NULL; /* LCOV_EXCL_LINE */
2862 : }
2863 :
2864 : /********************************************************************/
2865 : /** **/
2866 : /** SIMPLE MULTIPLICATION **/
2867 : /** **/
2868 : /********************************************************************/
2869 : GEN
2870 229578545 : gmulsg(long s, GEN x)
2871 : {
2872 : pari_sp av;
2873 : long i;
2874 : GEN z;
2875 :
2876 229578545 : switch(typ(x))
2877 : {
2878 136899439 : case t_INT: return mulsi(s,x);
2879 69367818 : case t_REAL: return s? mulsr(s,x): gen_0; /* gmul semantic */
2880 174122 : case t_INTMOD: { GEN p = gel(x,1);
2881 174122 : z = cgetg(3,t_INTMOD);
2882 174123 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(x,2)), p));
2883 174122 : gel(z,1) = icopy(p); return z;
2884 : }
2885 548217 : case t_FFELT: return FF_Z_mul(x,stoi(s));
2886 6207869 : case t_FRAC:
2887 6207869 : if (!s) return gen_0;
2888 6193225 : z = cgetg(3,t_FRAC);
2889 6193225 : i = labs(s); i = ugcd(i, umodiu(gel(x,2), i));
2890 6193225 : if (i == 1)
2891 : {
2892 2268166 : gel(z,2) = icopy(gel(x,2));
2893 2268166 : gel(z,1) = mulis(gel(x,1), s);
2894 : }
2895 : else
2896 : {
2897 3925059 : gel(z,2) = diviuexact(gel(x,2), (ulong)i);
2898 3925059 : gel(z,1) = mulis(gel(x,1), s/i);
2899 3925059 : fix_frac_if_int(z);
2900 : }
2901 6193225 : return z;
2902 :
2903 14202940 : case t_COMPLEX:
2904 14202940 : if (!s) return gen_0;
2905 13077136 : z = cgetg(3, t_COMPLEX);
2906 13077096 : gel(z,1) = gmulsg(s,gel(x,1));
2907 13075828 : gel(z,2) = gmulsg(s,gel(x,2)); return z;
2908 :
2909 1435 : case t_PADIC:
2910 1435 : if (!s) return gen_0;
2911 1435 : av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),x), x));
2912 :
2913 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2914 7 : gel(z,1) = ZX_copy(gel(x,1));
2915 7 : gel(z,2) = gmulsg(s,gel(x,2));
2916 7 : gel(z,3) = gmulsg(s,gel(x,3)); return z;
2917 :
2918 205408 : case t_POLMOD:
2919 205408 : retmkpolmod(gmulsg(s,gel(x,2)), RgX_copy(gel(x,1)));
2920 :
2921 746913 : case t_POL:
2922 746913 : if (!signe(x)) return RgX_copy(x);
2923 730582 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2924 3085461 : pari_APPLY_pol(gmulsg(s,gel(x,i)));
2925 :
2926 182 : case t_SER:
2927 182 : if (ser_isexactzero(x)) return gcopy(x);
2928 182 : if (!s) return Rg_get_0(x);
2929 3864 : pari_APPLY_ser(gmulsg(s,gel(x,i)));
2930 :
2931 0 : case t_RFRAC:
2932 0 : if (!s) return zeropol(varn(gel(x,2)));
2933 0 : if (s == 1) return gcopy(x);
2934 0 : if (s == -1) return gneg(x);
2935 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), stoi(s));
2936 :
2937 1233502 : case t_VEC: case t_COL: case t_MAT:
2938 3892935 : pari_APPLY_same(gmulsg(s,gel(x,i)));
2939 : }
2940 0 : pari_err_TYPE("gmulsg",x);
2941 : return NULL; /* LCOV_EXCL_LINE */
2942 : }
2943 :
2944 : GEN
2945 124338892 : gmulug(ulong s, GEN x)
2946 : {
2947 : pari_sp av;
2948 : long i;
2949 : GEN z;
2950 :
2951 124338892 : switch(typ(x))
2952 : {
2953 122223411 : case t_INT: return mului(s,x);
2954 1059371 : case t_REAL: return s? mulur(s,x): gen_0; /* gmul semantic */
2955 364 : case t_INTMOD: { GEN p = gel(x,1);
2956 364 : z = cgetg(3,t_INTMOD);
2957 364 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(mului(s,gel(x,2)), p));
2958 364 : gel(z,1) = icopy(p); return z;
2959 : }
2960 413 : case t_FFELT: return FF_Z_mul(x,utoi(s));
2961 979933 : case t_FRAC:
2962 979933 : if (!s) return gen_0;
2963 979919 : z = cgetg(3,t_FRAC);
2964 979919 : i = ugcd(s, umodiu(gel(x,2), s));
2965 979919 : if (i == 1)
2966 : {
2967 821635 : gel(z,2) = icopy(gel(x,2));
2968 821635 : gel(z,1) = muliu(gel(x,1), s);
2969 : }
2970 : else
2971 : {
2972 158284 : gel(z,2) = diviuexact(gel(x,2), i);
2973 158284 : gel(z,1) = muliu(gel(x,1), s/i);
2974 158284 : fix_frac_if_int(z);
2975 : }
2976 979919 : return z;
2977 :
2978 30765 : case t_COMPLEX:
2979 30765 : if (!s) return gen_0;
2980 30765 : z = cgetg(3, t_COMPLEX);
2981 30765 : gel(z,1) = gmulug(s,gel(x,1));
2982 30765 : gel(z,2) = gmulug(s,gel(x,2)); return z;
2983 :
2984 7342 : case t_PADIC:
2985 7342 : if (!s) return gen_0;
2986 7342 : av = avma; return gerepileupto(av, mulpp(cvtop2(utoi(s),x), x));
2987 :
2988 0 : case t_QUAD: z = cgetg(4, t_QUAD);
2989 0 : gel(z,1) = ZX_copy(gel(x,1));
2990 0 : gel(z,2) = gmulug(s,gel(x,2));
2991 0 : gel(z,3) = gmulug(s,gel(x,3)); return z;
2992 :
2993 6783 : case t_POLMOD:
2994 6783 : retmkpolmod(gmulug(s,gel(x,2)), RgX_copy(gel(x,1)));
2995 :
2996 16226 : case t_POL:
2997 16226 : if (!signe(x)) return RgX_copy(x);
2998 15449 : if (!s) return scalarpol(Rg_get_0(x), varn(x));
2999 52843 : pari_APPLY_pol(gmulug(s,gel(x,i)));
3000 :
3001 0 : case t_SER:
3002 0 : if (ser_isexactzero(x)) return gcopy(x);
3003 0 : if (!s) return Rg_get_0(x);
3004 0 : pari_APPLY_ser(gmulug(s,gel(x,i)));
3005 :
3006 0 : case t_RFRAC:
3007 0 : if (!s) return zeropol(varn(gel(x,2)));
3008 0 : if (s == 1) return gcopy(x);
3009 0 : return mul_rfrac_scal(gel(x,1), gel(x,2), utoi(s));
3010 :
3011 14285 : case t_VEC: case t_COL: case t_MAT:
3012 74464 : pari_APPLY_same(gmulug(s,gel(x,i)));
3013 : }
3014 0 : pari_err_TYPE("gmulsg",x);
3015 : return NULL; /* LCOV_EXCL_LINE */
3016 : }
3017 :
3018 : /********************************************************************/
3019 : /** **/
3020 : /** SIMPLE DIVISION **/
3021 : /** **/
3022 : /********************************************************************/
3023 :
3024 : GEN
3025 12932882 : gdivgs(GEN x, long s)
3026 : {
3027 12932882 : long tx = typ(x), i;
3028 : pari_sp av;
3029 : GEN z;
3030 :
3031 12932882 : if (!s)
3032 : {
3033 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3034 0 : pari_err_INV("gdivgs",gen_0);
3035 : }
3036 12932888 : switch(tx)
3037 : {
3038 1560427 : case t_INT: return Qdivis(x, s);
3039 8495323 : case t_REAL: return divrs(x,s);
3040 :
3041 357 : case t_INTMOD:
3042 357 : z = cgetg(3, t_INTMOD);
3043 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3044 :
3045 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3046 :
3047 552128 : case t_FRAC: z = cgetg(3, t_FRAC);
3048 552128 : i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
3049 552128 : if (i == 1)
3050 : {
3051 409305 : gel(z,2) = mulsi(s, gel(x,2));
3052 409305 : gel(z,1) = icopy(gel(x,1));
3053 : }
3054 : else
3055 : {
3056 142823 : gel(z,2) = mulsi(s/i, gel(x,2));
3057 142823 : gel(z,1) = divis(gel(x,1), i);
3058 : }
3059 552128 : normalize_frac(z);
3060 552128 : fix_frac_if_int(z); return z;
3061 :
3062 1840625 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3063 1840624 : gel(z,1) = gdivgs(gel(x,1),s);
3064 1840621 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3065 :
3066 133 : case t_PADIC: /* divpT */
3067 : {
3068 133 : GEN p = gel(x,2);
3069 133 : if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
3070 133 : av = avma;
3071 133 : return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
3072 : }
3073 :
3074 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3075 28 : gel(z,1) = ZX_copy(gel(x,1));
3076 28 : gel(z,2) = gdivgs(gel(x,2),s);
3077 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3078 :
3079 37826 : case t_POLMOD:
3080 37826 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3081 :
3082 91 : case t_RFRAC:
3083 91 : if (s == 1) return gcopy(x);
3084 84 : else if (s == -1) return gneg(x);
3085 84 : return div_rfrac_scal(x, stoi(s));
3086 :
3087 157457 : case t_POL: pari_APPLY_pol_normalized(gdivgs(gel(x,i),s));
3088 0 : case t_SER: pari_APPLY_ser_normalized(gdivgs(gel(x,i),s));
3089 407593 : case t_VEC:
3090 : case t_COL:
3091 894977 : case t_MAT: pari_APPLY_same(gdivgs(gel(x,i),s));
3092 : }
3093 0 : pari_err_TYPE2("/",x, stoi(s));
3094 : return NULL; /* LCOV_EXCL_LINE */
3095 : }
3096 :
3097 : GEN
3098 54712757 : gdivgu(GEN x, ulong s)
3099 : {
3100 54712757 : long tx = typ(x), i;
3101 : pari_sp av;
3102 : GEN z;
3103 :
3104 54712757 : if (!s)
3105 : {
3106 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3107 0 : pari_err_INV("gdivgu",gen_0);
3108 : }
3109 54712684 : switch(tx)
3110 : {
3111 17915794 : case t_INT: return Qdiviu(x, s);
3112 12954300 : case t_REAL: return divru(x,s);
3113 :
3114 210315 : case t_INTMOD:
3115 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3116 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3117 :
3118 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3119 :
3120 641982 : case t_FRAC: z = cgetg(3, t_FRAC);
3121 641982 : i = ugcd(s, umodiu(gel(x,1), s));
3122 641982 : if (i == 1)
3123 : {
3124 462120 : gel(z,2) = mului(s, gel(x,2));
3125 462120 : gel(z,1) = icopy(gel(x,1));
3126 : }
3127 : else
3128 : {
3129 179862 : gel(z,2) = mului(s/i, gel(x,2));
3130 179862 : gel(z,1) = divis(gel(x,1), i);
3131 : }
3132 641982 : normalize_frac(z);
3133 641982 : fix_frac_if_int(z); return z;
3134 :
3135 11383826 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3136 11383827 : gel(z,1) = gdivgu(gel(x,1),s);
3137 11383826 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3138 :
3139 11551454 : case t_PADIC: /* divpT */
3140 : {
3141 11551454 : GEN p = gel(x,2);
3142 11551454 : if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
3143 11415994 : av = avma;
3144 11415994 : return gerepileupto(av, divpp(x, cvtop2(utoi(s),x)));
3145 : }
3146 :
3147 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3148 0 : gel(z,1) = ZX_copy(gel(x,1));
3149 0 : gel(z,2) = gdivgu(gel(x,2),s);
3150 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3151 :
3152 1456 : case t_POLMOD:
3153 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3154 :
3155 56 : case t_RFRAC:
3156 56 : if (s == 1) return gcopy(x);
3157 56 : return div_rfrac_scal(x, utoi(s));
3158 :
3159 205660 : case t_POL: pari_APPLY_pol_normalized(gdivgu(gel(x,i),s));
3160 16730 : case t_SER: pari_APPLY_ser_normalized(gdivgu(gel(x,i),s));
3161 329 : case t_VEC:
3162 : case t_COL:
3163 1148 : case t_MAT: pari_APPLY_same(gdivgu(gel(x,i),s));
3164 : }
3165 0 : pari_err_TYPE2("/",x, utoi(s));
3166 : return NULL; /* LCOV_EXCL_LINE */
3167 : }
3168 :
3169 : /* x / (i*(i+1)) */
3170 : GEN
3171 224165928 : divrunextu(GEN x, ulong i)
3172 : {
3173 224165928 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3174 0 : return divri(x, muluu(i , i+1));
3175 : else
3176 224165928 : return divru(x, i*(i+1));
3177 : }
3178 : /* x / (i*(i+1)) */
3179 : GEN
3180 814832 : gdivgunextu(GEN x, ulong i)
3181 : {
3182 814832 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3183 0 : return gdivgu(x, i*(i+1));
3184 : else
3185 814832 : return gdiv(x, muluu(i, i+1));
3186 : }
3187 :
3188 : /* True shift (exact multiplication by 2^n) */
3189 : GEN
3190 124671379 : gmul2n(GEN x, long n)
3191 : {
3192 : GEN z, a, b;
3193 : long k, l;
3194 :
3195 124671379 : switch(typ(x))
3196 : {
3197 31203587 : case t_INT:
3198 31203587 : if (n>=0) return shifti(x,n);
3199 5035100 : if (!signe(x)) return gen_0;
3200 3341651 : l = vali(x); n = -n;
3201 3341711 : if (n<=l) return shifti(x,-n);
3202 398367 : z = cgetg(3,t_FRAC);
3203 398367 : gel(z,1) = shifti(x,-l);
3204 398367 : gel(z,2) = int2n(n-l); return z;
3205 :
3206 62790992 : case t_REAL:
3207 62790992 : return shiftr(x,n);
3208 :
3209 180971 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3210 180971 : z = cgetg(3,t_INTMOD);
3211 180971 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3212 82810 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
3213 82809 : gel(z,1) = icopy(b); return z;
3214 :
3215 217417 : case t_FFELT: return FF_mul2n(x,n);
3216 :
3217 1268257 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3218 1268257 : l = vali(a);
3219 1268257 : k = vali(b);
3220 1268257 : if (n+l >= k)
3221 : {
3222 395200 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3223 302690 : l = n-k; k = -k;
3224 : }
3225 : else
3226 : {
3227 873057 : k = -(l+n); l = -l;
3228 : }
3229 1175747 : z = cgetg(3,t_FRAC);
3230 1175747 : gel(z,1) = shifti(a,l);
3231 1175747 : gel(z,2) = shifti(b,k); return z;
3232 :
3233 10726793 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3234 10726792 : gel(z,1) = gmul2n(gel(x,1),n);
3235 10726792 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3236 :
3237 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3238 105 : gel(z,1) = ZX_copy(gel(x,1));
3239 105 : gel(z,2) = gmul2n(gel(x,2),n);
3240 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3241 :
3242 193970 : case t_POLMOD:
3243 193970 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3244 :
3245 1717798 : case t_POL:
3246 9105332 : pari_APPLY_pol(gmul2n(gel(x,i),n));
3247 102847 : case t_SER:
3248 102847 : if (ser_isexactzero(x)) return gcopy(x);
3249 634256 : pari_APPLY_ser(gmul2n(gel(x,i),n));
3250 16265058 : case t_VEC: case t_COL: case t_MAT:
3251 59926017 : pari_APPLY_same(gmul2n(gel(x,i),n));
3252 :
3253 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3254 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3255 :
3256 5579 : case t_PADIC: /* int2n wrong if n < 0 */
3257 5579 : return gmul(gmul2n(gen_1,n),x);
3258 : }
3259 0 : pari_err_TYPE("gmul2n",x);
3260 : return NULL; /* LCOV_EXCL_LINE */
3261 : }
3262 :
3263 : /*******************************************************************/
3264 : /* */
3265 : /* INVERSE */
3266 : /* */
3267 : /*******************************************************************/
3268 : static GEN
3269 209318 : inv_polmod(GEN T, GEN x)
3270 : {
3271 209318 : GEN z = cgetg(3,t_POLMOD), a;
3272 209318 : gel(z,1) = RgX_copy(T);
3273 209317 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3274 84684 : a = ginv(x);
3275 : else
3276 : {
3277 124633 : if (lg(T) == 5) /* quadratic fields */
3278 13055 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3279 : else
3280 111578 : a = RgXQ_inv(x, T);
3281 : }
3282 209317 : gel(z,2) = a; return z;
3283 : }
3284 :
3285 : GEN
3286 36431655 : ginv(GEN x)
3287 : {
3288 : long s;
3289 : pari_sp av, tetpil;
3290 : GEN z, y, p1, p2;
3291 :
3292 36431655 : switch(typ(x))
3293 : {
3294 9450835 : case t_INT:
3295 9450835 : if (is_pm1(x)) return icopy(x);
3296 7896002 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3297 7895988 : z = cgetg(3,t_FRAC);
3298 7900858 : gel(z,1) = s<0? gen_m1: gen_1;
3299 7900858 : gel(z,2) = absi(x); return z;
3300 :
3301 10741909 : case t_REAL: return invr(x);
3302 :
3303 12964 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3304 12964 : gel(z,1) = icopy(gel(x,1));
3305 12964 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3306 :
3307 552058 : case t_FRAC: {
3308 552058 : GEN a = gel(x,1), b = gel(x,2);
3309 552058 : s = signe(a);
3310 552058 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3311 248264 : z = cgetg(3,t_FRAC);
3312 248264 : gel(z,1) = icopy(b);
3313 248263 : gel(z,2) = icopy(a);
3314 248264 : normalize_frac(z); return z;
3315 : }
3316 9601781 : case t_COMPLEX:
3317 9601781 : av=avma;
3318 9601781 : p1=cxnorm(x);
3319 9600888 : p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
3320 9601578 : tetpil=avma;
3321 9601578 : return gerepile(av,tetpil,divcR(p2,p1));
3322 :
3323 273 : case t_QUAD:
3324 273 : av=avma; p1=quadnorm(x); p2=conj_i(x); tetpil=avma;
3325 273 : return gerepile(av,tetpil,gdiv(p2,p1));
3326 :
3327 358358 : case t_PADIC: z = cgetg(5,t_PADIC);
3328 358358 : if (!signe(gel(x,4))) pari_err_INV("ginv",x);
3329 358351 : z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
3330 358351 : gel(z,2) = icopy(gel(x,2));
3331 358351 : gel(z,3) = icopy(gel(x,3));
3332 358351 : gel(z,4) = Zp_inv(gel(x,4),gel(z,2),precp(x)); return z;
3333 :
3334 209318 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3335 14221 : case t_FFELT: return FF_inv(x);
3336 5322842 : case t_POL: return gred_rfrac_simple(gen_1,x);
3337 26341 : case t_SER: return ser_inv(x);
3338 2933 : case t_RFRAC:
3339 : {
3340 2933 : GEN n = gel(x,1), d = gel(x,2);
3341 2933 : pari_sp av = avma, ltop;
3342 2933 : if (gequal0(n)) pari_err_INV("ginv",x);
3343 :
3344 2933 : n = simplify_shallow(n);
3345 2933 : if (typ(n) != t_POL || varn(n) != varn(d))
3346 : {
3347 2933 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3348 679 : ltop = avma;
3349 679 : z = RgX_Rg_div(d,n);
3350 : } else {
3351 0 : ltop = avma;
3352 0 : z = cgetg(3,t_RFRAC);
3353 0 : gel(z,1) = RgX_copy(d);
3354 0 : gel(z,2) = RgX_copy(n);
3355 : }
3356 679 : stackdummy(av, ltop);
3357 679 : return z;
3358 : }
3359 :
3360 21 : case t_VEC: if (!is_ext_qfr(x)) break;
3361 : case t_QFB:
3362 99599 : return qfbpow(x, gen_m1);
3363 38820 : case t_MAT:
3364 38820 : y = RgM_inv(x);
3365 38813 : if (!y) pari_err_INV("ginv",x);
3366 38743 : return y;
3367 28 : case t_VECSMALL:
3368 : {
3369 28 : long i, lx = lg(x)-1;
3370 28 : y = zero_zv(lx);
3371 112 : for (i=1; i<=lx; i++)
3372 : {
3373 84 : long xi = x[i];
3374 84 : if (xi<1 || xi>lx || y[xi])
3375 0 : pari_err_TYPE("ginv [not a permutation]", x);
3376 84 : y[xi] = i;
3377 : }
3378 28 : return y;
3379 : }
3380 : }
3381 6 : pari_err_TYPE("inverse",x);
3382 : return NULL; /* LCOV_EXCL_LINE */
3383 : }
|
