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 98 : warn_coercion(GEN x, GEN y, GEN z)
38 : {
39 98 : if (DEBUGLEVEL)
40 56 : pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
41 98 : }
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 216795396 : addRc(GEN x, GEN y) {
61 216795396 : GEN z = cgetg(3,t_COMPLEX);
62 216794613 : gel(z,1) = gadd(x,gel(y,1));
63 216740524 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 266442132 : mulRc(GEN x, GEN y) {
67 266442132 : GEN z = cgetg(3,t_COMPLEX);
68 266447362 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
69 266415934 : 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 2512625 : divRc(GEN x, GEN y) {
80 2512625 : GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
81 2512624 : GEN z = cgetg(3,t_COMPLEX);
82 2512624 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
83 2512600 : gel(z,2) = gmul(mt, gel(y,2));
84 2512601 : return z;
85 : }
86 : static GEN
87 18609208 : divcR(GEN x, GEN y) {
88 18609208 : GEN z = cgetg(3,t_COMPLEX);
89 18609124 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
90 18608188 : gel(z,2) = gdiv(gel(x,2), y); return z;
91 : }
92 : static GEN
93 1267 : addRq(GEN x, GEN y) {
94 1267 : GEN z = cgetg(4,t_QUAD);
95 1267 : gel(z,1) = ZX_copy(gel(y,1));
96 1267 : gel(z,2) = gadd(x, gel(y,2));
97 1267 : 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 70 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
108 70 : long i = gexpo(x) - gexpo(y);
109 70 : if (i > 0) prec += nbits2extraprec(i);
110 70 : return gerepileupto(av, gadd(y, quadtofp(x, prec)));
111 : }
112 : static GEN
113 29262851 : mulrfrac(GEN x, GEN y)
114 : {
115 : pari_sp av;
116 29262851 : GEN z, a = gel(y,1), b = gel(y,2);
117 29262851 : if (is_pm1(a)) /* frequent special case */
118 : {
119 5242059 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
120 5241913 : return z;
121 : }
122 24020800 : av = avma;
123 24020800 : 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 6496826 : mulTp(GEN x, GEN y) { pari_sp av = avma;
151 6496826 : 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 3647 : divTp(GEN x, GEN y) { pari_sp av = avma;
156 3647 : 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 1327200 : divpT(GEN x, GEN y) { pari_sp av = avma;
162 1327200 : 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 1974870 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
169 1974870 : if (lgefint(X) == 3) {
170 1894040 : ulong u = Fl_add(itou(x),itou(y), X[2]);
171 1894040 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
172 : }
173 : else {
174 80830 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
175 80829 : gel(z,2) = gerepileuptoint((pari_sp)z, u);
176 : }
177 1974871 : gel(z,1) = icopy(X); return z;
178 : }
179 : static GEN
180 1153254 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 1153254 : if (lgefint(X) == 3) {
182 779220 : ulong u = Fl_sub(itou(x),itou(y), X[2]);
183 779220 : 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 1153254 : gel(z,1) = icopy(X); return z;
190 : }
191 : /* cf add_intmod_same */
192 : static GEN
193 3218325 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
194 3218325 : if (lgefint(X) == 3) {
195 3109932 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
196 3109932 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
197 : }
198 : else
199 108393 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
200 3218334 : gel(z,1) = icopy(X); return z;
201 : }
202 : /* cf add_intmod_same */
203 : static GEN
204 337532 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
205 : {
206 337532 : if (lgefint(X) == 3) {
207 314719 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
208 314712 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 22813 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
212 337526 : 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 9884933 : rfrac_denom_mul_scal(GEN d, GEN y)
225 : {
226 9884933 : GEN D = RgX_Rg_mul(d, y);
227 9884933 : 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 9884933 : return D;
233 : }
234 :
235 : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
236 : GEN
237 57837375 : gred_rfrac_simple(GEN n, GEN d)
238 : {
239 : GEN c, cn, cd, z;
240 57837375 : long dd = degpol(d);
241 :
242 57837375 : if (dd <= 0)
243 : {
244 4123 : if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
245 4123 : n = gdiv(n, gel(d,2));
246 4123 : if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
247 4123 : return n;
248 : }
249 :
250 57833252 : cd = content(d);
251 59699419 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
252 57833252 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
253 57833252 : if (!gequal1(cd)) {
254 6582165 : d = RgX_Rg_div(d,cd);
255 6582165 : if (!gequal1(cn))
256 : {
257 1309528 : if (gequal0(cn)) {
258 49 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
259 0 : n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
260 0 : c = gen_1;
261 : } else {
262 1309479 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
263 1309479 : c = gdiv(cn,cd);
264 : }
265 : }
266 : else
267 5272637 : c = ginv(cd);
268 : } else {
269 51251087 : if (!gequal1(cn))
270 : {
271 3272895 : if (gequal0(cn)) {
272 931 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
273 21 : c = gen_1;
274 : } else {
275 3271964 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
276 3271964 : c = cn;
277 : }
278 : } else {
279 47978192 : GEN y = cgetg(3,t_RFRAC);
280 47978192 : gel(y,1) = gcopy(n);
281 47978192 : gel(y,2) = RgX_copy(d); return y;
282 : }
283 : }
284 :
285 9854101 : if (typ(c) == t_POL)
286 : {
287 904069 : z = c;
288 943213 : do { z = content(z); } while (typ(z) == t_POL);
289 904069 : cd = denom_i(z);
290 904069 : cn = gmul(c, cd);
291 : }
292 : else
293 : {
294 8950032 : cn = numer_i(c);
295 8950032 : cd = denom_i(c);
296 : }
297 9854101 : z = cgetg(3,t_RFRAC);
298 9854101 : gel(z,1) = gmul(n, cn);
299 9854101 : gel(z,2) = rfrac_denom_mul_scal(d, cd);
300 9854101 : return z;
301 : }
302 :
303 : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
304 : static GEN
305 249123 : fix_rfrac(GEN x, long d)
306 : {
307 : GEN z, N, D;
308 249123 : if (!d || typ(x) == t_POL) return x;
309 175511 : z = cgetg(3, t_RFRAC);
310 175511 : N = gel(x,1);
311 175511 : D = gel(x,2);
312 175511 : if (d > 0) {
313 11690 : gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
314 177842 : : monomialcopy(N,d,varn(D));
315 166152 : gel(z, 2) = RgX_copy(D);
316 : } else {
317 9359 : gel(z, 1) = gcopy(N);
318 9359 : gel(z, 2) = RgX_shift(D, -d);
319 : }
320 175511 : return z;
321 : }
322 :
323 : /* assume d != 0 */
324 : static GEN
325 44767912 : gred_rfrac2(GEN n, GEN d)
326 : {
327 : GEN y, z;
328 : long v, vd, vn;
329 :
330 44767912 : n = simplify_shallow(n);
331 44767912 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
332 37936360 : d = simplify_shallow(d);
333 37936360 : if (typ(d) != t_POL) return gdiv(n,d);
334 36756564 : vd = varn(d);
335 36756564 : if (typ(n) != t_POL)
336 : {
337 20622506 : if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
338 20621099 : if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
339 0 : pari_err_BUG("gred_rfrac2 [incompatible variables]");
340 : }
341 16134058 : vn = varn(n);
342 16134058 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
343 15996412 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
344 :
345 : /* now n and d are t_POLs in the same variable */
346 15836660 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
347 15836660 : if (!degpol(d))
348 : {
349 12760350 : n = RgX_Rg_div(n,gel(d,2));
350 12760350 : return v? RgX_mulXn(n,v): n;
351 : }
352 :
353 : /* X does not divide gcd(n,d), deg(d) > 0 */
354 3076310 : if (!isinexact(n) && !isinexact(d))
355 : {
356 3076093 : y = RgX_divrem(n, d, &z);
357 3076093 : if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
358 248906 : z = RgX_gcd(d, z);
359 248906 : if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
360 : }
361 249123 : return fix_rfrac(gred_rfrac_simple(n,d), v);
362 : }
363 :
364 : /* x,y t_INT, return x/y in reduced form */
365 : GEN
366 109262585 : Qdivii(GEN x, GEN y)
367 : {
368 109262585 : pari_sp av = avma;
369 : GEN r, q;
370 :
371 109262585 : if (lgefint(y) == 3)
372 : {
373 94169777 : q = Qdiviu(x, y[2]);
374 94166229 : if (signe(y) > 0) return q;
375 10672004 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
376 10672184 : return q;
377 : }
378 15092808 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
379 15094103 : if (equali1(x))
380 : {
381 5177462 : if (!signe(y)) pari_err_INV("gdiv",y);
382 5177350 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
383 : }
384 9916649 : q = dvmdii(x,y,&r);
385 9916402 : if (r == gen_0) return q; /* gen_0 intended */
386 7012052 : r = gcdii(y, r);
387 7011885 : if (lgefint(r) == 3)
388 : {
389 6184303 : ulong t = r[2];
390 6184303 : set_avma(av);
391 6184388 : if (t == 1) q = mkfraccopy(x,y);
392 : else
393 : {
394 2463879 : q = cgetg(3,t_FRAC);
395 2464262 : gel(q,1) = diviuexact(x,t);
396 2463957 : gel(q,2) = diviuexact(y,t);
397 : }
398 : }
399 : else
400 : { /* rare: r and q left on stack for efficiency */
401 827582 : q = cgetg(3,t_FRAC);
402 827620 : gel(q,1) = diviiexact(x,r);
403 827626 : gel(q,2) = diviiexact(y,r);
404 : }
405 7011813 : normalize_frac(q); return q;
406 : }
407 :
408 : /* x t_INT, return x/y in reduced form */
409 : GEN
410 112342243 : Qdiviu(GEN x, ulong y)
411 : {
412 112342243 : pari_sp av = avma;
413 : ulong r, t;
414 : GEN q;
415 :
416 112342243 : if (y == 1) return icopy(x);
417 93151721 : if (!y) pari_err_INV("gdiv",gen_0);
418 93151721 : if (equali1(x)) retmkfrac(gen_1, utoipos(y));
419 87468311 : q = absdiviu_rem(x,y,&r);
420 87460479 : if (!r)
421 : {
422 48985498 : if (signe(x) < 0) togglesign(q);
423 48986523 : return q;
424 : }
425 38474981 : t = ugcd(y, r); set_avma(av);
426 38481479 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
427 15758744 : retmkfrac(diviuexact(x,t), utoipos(y / t));
428 : }
429 :
430 : /* x t_INT, return x/y in reduced form */
431 : GEN
432 1326121 : Qdivis(GEN x, long y)
433 : {
434 1326121 : pari_sp av = avma;
435 : ulong r, t;
436 : long s;
437 : GEN q;
438 :
439 1326121 : if (y > 0) return Qdiviu(x, y);
440 98434 : if (!y) pari_err_INV("gdiv",gen_0);
441 98434 : s = signe(x);
442 98434 : if (!s) return gen_0;
443 70140 : if (y < 0) { y = -y; s = -s; }
444 70140 : if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
445 69860 : if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
446 67830 : q = absdiviu_rem(x,y,&r);
447 67830 : if (!r)
448 : {
449 38402 : if (s < 0) togglesign(q);
450 38402 : return q;
451 : }
452 29428 : t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
453 29428 : if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
454 29428 : gel(q,1) = x; setsigne(x, s);
455 29428 : gel(q,2) = utoipos(y); return q;
456 : }
457 :
458 : /*******************************************************************/
459 : /* */
460 : /* CONJUGATION */
461 : /* */
462 : /*******************************************************************/
463 : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
464 : static GEN
465 17850 : quad_polmod_conj(GEN x, GEN y)
466 : {
467 : GEN z, u, v, a, b;
468 17850 : if (typ(x) != t_POL) return gcopy(x);
469 17850 : if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
470 17850 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
471 17850 : b = gel(y,3); v = gel(x,2);
472 17850 : z = cgetg(4, t_POL); z[1] = x[1];
473 17850 : gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
474 17850 : gel(z,3) = gneg(u); return z;
475 : }
476 : static GEN
477 17850 : quad_polmod_norm(GEN x, GEN y)
478 : {
479 : GEN z, u, v, a, b, c;
480 17850 : if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
481 17850 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
482 17850 : b = gel(y,3); v = gel(x,2);
483 17850 : c = gel(y,2);
484 17850 : z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
485 17850 : if (!gequal1(a)) z = gdiv(z, a);
486 17850 : return gadd(z, gsqr(v));
487 : }
488 :
489 : GEN
490 17163436 : conj_i(GEN x)
491 : {
492 : long lx, i;
493 : GEN y;
494 :
495 17163436 : switch(typ(x))
496 : {
497 5783852 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
498 5783852 : return x;
499 :
500 11220800 : case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
501 :
502 924 : case t_QUAD:
503 924 : y = cgetg(4,t_QUAD);
504 924 : gel(y,1) = gel(x,1);
505 924 : gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
506 924 : : gadd(gel(x,2), gel(x,3));
507 924 : gel(y,3) = gneg(gel(x,3)); return y;
508 :
509 8708 : case t_POL: case t_SER:
510 8708 : y = cgetg_copy(x, &lx); y[1] = x[1];
511 29323 : for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
512 8708 : return y;
513 :
514 149184 : case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
515 149184 : y = cgetg_copy(x, &lx);
516 538094 : for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
517 149183 : return y;
518 :
519 0 : case t_POLMOD:
520 : {
521 0 : GEN X = gel(x,1);
522 0 : long d = degpol(X);
523 0 : if (d < 2) return x;
524 0 : if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
525 : }
526 : }
527 0 : pari_err_TYPE("gconj",x);
528 : return NULL; /* LCOV_EXCL_LINE */
529 : }
530 : GEN
531 95388 : gconj(GEN x)
532 95388 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
533 :
534 : GEN
535 84 : conjvec(GEN x,long prec)
536 : {
537 : long lx, s, i;
538 : GEN z;
539 :
540 84 : switch(typ(x))
541 : {
542 0 : case t_INT: case t_INTMOD: case t_FRAC:
543 0 : return mkcolcopy(x);
544 :
545 0 : case t_COMPLEX: case t_QUAD:
546 0 : z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
547 :
548 28 : case t_FFELT:
549 28 : return FF_conjvec(x);
550 :
551 0 : case t_VEC: case t_COL:
552 0 : lx = lg(x); z = cgetg(lx,t_MAT);
553 0 : if (lx == 1) return z;
554 0 : gel(z,1) = conjvec(gel(x,1),prec);
555 0 : s = lgcols(z);
556 0 : for (i=2; i<lx; i++)
557 : {
558 0 : gel(z,i) = conjvec(gel(x,i),prec);
559 0 : if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
560 : }
561 0 : break;
562 :
563 56 : case t_POLMOD: {
564 56 : GEN T = gel(x,1), r;
565 : pari_sp av;
566 :
567 56 : lx = lg(T);
568 56 : if (lx <= 3) return cgetg(1,t_COL);
569 56 : x = gel(x,2);
570 238 : for (i=2; i<lx; i++)
571 : {
572 189 : GEN c = gel(T,i);
573 189 : switch(typ(c)) {
574 7 : case t_INTMOD: {
575 7 : GEN p = gel(c,1);
576 : pari_sp av;
577 7 : if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
578 7 : av = avma;
579 7 : T = RgX_to_FpX(T,p);
580 7 : x = RgX_to_FpX(x, p);
581 7 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
582 7 : z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
583 7 : return gerepileupto(av, z);
584 : }
585 182 : case t_INT:
586 182 : case t_FRAC: break;
587 0 : default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
588 : }
589 : }
590 49 : if (typ(x) != t_POL)
591 : {
592 0 : if (!is_rational_t(typ(x)))
593 0 : pari_err_TYPE("conjvec [not a rational t_POL]",x);
594 0 : retconst_col(lx-3, gcopy(x));
595 : }
596 49 : RgX_check_QX(x,"conjvec");
597 49 : av = avma;
598 49 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
599 49 : r = cleanroots(T,prec);
600 49 : z = cgetg(lx-2,t_COL);
601 182 : for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
602 49 : return gerepileupto(av, z);
603 : }
604 :
605 0 : default:
606 0 : pari_err_TYPE("conjvec",x);
607 : return NULL; /* LCOV_EXCL_LINE */
608 : }
609 0 : return z;
610 : }
611 :
612 : /********************************************************************/
613 : /** **/
614 : /** ADDITION **/
615 : /** **/
616 : /********************************************************************/
617 : /* x, y compatible PADIC, op = add or sub */
618 : static GEN
619 10531647 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
620 : {
621 10531647 : pari_sp av = avma;
622 : long d,e,r,rx,ry;
623 10531647 : GEN u, z, mod, p = gel(x,2);
624 : int swap;
625 :
626 10531647 : (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
627 10531385 : e = valp(x);
628 10531385 : r = valp(y); d = r-e;
629 10531385 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
630 10531385 : rx = precp(x);
631 10531385 : ry = precp(y);
632 10531385 : if (d) /* v(x) < v(y) */
633 : {
634 7368739 : r = d+ry; z = powiu(p,d);
635 7368823 : if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
636 7368840 : z = mulii(z,gel(y,4));
637 7368787 : u = swap? op(z, gel(x,4)): op(gel(x,4), z);
638 : }
639 : else
640 : {
641 : long c;
642 3162646 : if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
643 3162646 : u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
644 3163446 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
645 : {
646 138619 : set_avma(av); return zeropadic(p, e+r);
647 : }
648 3025098 : if (c)
649 : {
650 1870118 : mod = diviiexact(mod, powiu(p,c));
651 1870115 : r -= c;
652 1870115 : e += c;
653 : }
654 : }
655 10393850 : u = modii(u, mod);
656 10392873 : set_avma(av); z = cgetg(5,t_PADIC);
657 10392210 : z[1] = evalprecp(r) | evalvalp(e);
658 10392160 : gel(z,2) = icopy(p);
659 10391998 : gel(z,3) = icopy(mod);
660 10392057 : gel(z,4) = icopy(u); return z;
661 : }
662 : /* Rg_to_Fp(t_FRAC) without GC */
663 : static GEN
664 23555 : Q_to_Fp(GEN x, GEN p)
665 23555 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
666 : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
667 : static GEN
668 3513464 : addQp(GEN x, GEN y)
669 : {
670 3513464 : pari_sp av = avma;
671 3513464 : long d, r, e, vy = valp(y), py = precp(y);
672 3513464 : GEN z, q, mod, u, p = gel(y,2);
673 :
674 3513464 : e = Q_pvalrem(x, p, &x);
675 3513450 : d = vy - e; r = d + py;
676 3513450 : if (r <= 0) { set_avma(av); return gcopy(y); }
677 3511497 : mod = gel(y,3);
678 3511497 : u = gel(y,4);
679 3511497 : (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
680 :
681 3511491 : if (d > 0)
682 : {
683 893989 : q = powiu(p,d);
684 894075 : mod = mulii(mod, q);
685 894075 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
686 894075 : u = addii(x, mulii(u, q));
687 : }
688 2617502 : else if (d < 0)
689 : {
690 397578 : q = powiu(p,-d);
691 397578 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
692 397578 : u = addii(u, mulii(x, q));
693 397578 : r = py; e = vy;
694 : }
695 : else
696 : {
697 : long c;
698 2219924 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
699 2219924 : u = addii(u, x);
700 2219909 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
701 : {
702 904 : set_avma(av); return zeropadic(p,e+r);
703 : }
704 2219012 : if (c)
705 : {
706 796797 : mod = diviiexact(mod, powiu(p,c));
707 796797 : r -= c;
708 796797 : e += c;
709 : }
710 : }
711 3510668 : u = modii(u, mod); set_avma(av);
712 3510609 : z = cgetg(5,t_PADIC);
713 3510560 : z[1] = evalprecp(r) | evalvalp(e);
714 3510555 : gel(z,2) = icopy(p);
715 3510536 : gel(z,3) = icopy(mod);
716 3510534 : gel(z,4) = icopy(u); return z;
717 : }
718 :
719 : /* Mod(x,X) + Mod(y,X) */
720 : #define addsub_polmod_same addsub_polmod_scal
721 : /* Mod(x,X) +/- Mod(y,Y) */
722 : static GEN
723 2142 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
724 : {
725 2142 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
726 2142 : GEN z = cgetg(3,t_POLMOD);
727 2142 : long vx = varn(X), vy = varn(Y);
728 2142 : if (vx==vy) {
729 : pari_sp av;
730 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
731 14 : warn_coercion(X,Y,gel(z,1));
732 14 : gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
733 : }
734 2128 : if (varncmp(vx, vy) < 0)
735 2128 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
736 : else
737 0 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
738 2128 : gel(z,2) = op(x, y); return z;
739 : }
740 : /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
741 : static GEN
742 12527896 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
743 12527896 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
744 :
745 : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
746 : static GEN
747 369328 : add_ser_scal(GEN y, GEN x)
748 : {
749 : long i, v, ly, vy;
750 : GEN z;
751 :
752 369328 : if (isrationalzero(x)) return gcopy(y);
753 344275 : ly = lg(y);
754 344275 : v = valp(y);
755 344275 : if (v < 3-ly) return gcopy(y);
756 : /* v + ly >= 3 */
757 344051 : if (v < 0)
758 : {
759 1169 : z = cgetg(ly,t_SER); z[1] = y[1];
760 3304 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
761 1169 : gel(z,i) = gadd(x,gel(y,i)); i++;
762 3276 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
763 1169 : return normalizeser(z);
764 : }
765 342882 : vy = varn(y);
766 342882 : if (v > 0)
767 : {
768 15946 : if (ser_isexactzero(y))
769 7406 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
770 8540 : y -= v; ly += v;
771 8540 : z = cgetg(ly,t_SER);
772 8540 : x = gcopy(x);
773 19138 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
774 : }
775 326936 : else if (ser_isexactzero(y)) return gcopy(y);
776 : else
777 : { /* v = 0, ly >= 3 */
778 326929 : z = cgetg(ly,t_SER);
779 326929 : x = gadd(x, gel(y,2));
780 326929 : i = 3;
781 : }
782 1508196 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
783 335469 : gel(z,2) = x;
784 335469 : z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
785 335469 : return gequal0(x)? normalizeser(z): z;
786 : }
787 : static long
788 6921726 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
789 : /* x,y t_SER in the same variable: x+y */
790 : static GEN
791 3461262 : ser_add(GEN x, GEN y)
792 : {
793 3461262 : long i, lx,ly, n = valp(y) - valp(x);
794 : GEN z;
795 3461262 : if (n < 0) { n = -n; swap(x,y); }
796 : /* valp(x) <= valp(y) */
797 3461262 : lx = _serprec(x);
798 3461262 : if (lx == 2) /* don't lose type information */
799 : {
800 798 : z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
801 798 : setvalp(z, valp(x)); return z;
802 : }
803 3460464 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
804 3460464 : if (n)
805 : {
806 114887 : if (n+2 > lx) return gcopy(x);
807 114313 : z = cgetg(ly,t_SER);
808 841684 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
809 529248 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
810 : } else {
811 3345577 : z = cgetg(ly,t_SER);
812 16837374 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
813 : }
814 3459890 : z[1] = x[1]; return normalizeser(z);
815 : }
816 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
817 : static GEN
818 8824491 : add_rfrac_scal(GEN y, GEN x)
819 : {
820 : pari_sp av;
821 : GEN n;
822 :
823 8824491 : if (isintzero(x)) return gcopy(y); /* frequent special case */
824 5150726 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
825 5150726 : return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
826 : }
827 :
828 : /* x "scalar", ty != t_MAT and nonscalar */
829 : static GEN
830 21229791 : add_scal(GEN y, GEN x, long ty)
831 : {
832 21229791 : switch(ty)
833 : {
834 16377606 : case t_POL: return RgX_Rg_add(y, x);
835 369307 : case t_SER: return add_ser_scal(y, x);
836 4468757 : case t_RFRAC: return add_rfrac_scal(y, x);
837 0 : case t_COL: return RgC_Rg_add(y, x);
838 14120 : case t_VEC:
839 14120 : if (isintzero(x)) return gcopy(y);
840 182 : break;
841 : }
842 183 : pari_err_TYPE2("+",x,y);
843 : return NULL; /* LCOV_EXCL_LINE */
844 : }
845 :
846 : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
847 : static GEN
848 12571103 : setfrac(GEN z, GEN a, GEN b)
849 : {
850 12571103 : gel(z,1) = icopy_avma(a, (pari_sp)z);
851 12571131 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
852 12571209 : set_avma((pari_sp)gel(z,2)); return z;
853 : }
854 : /* z <- a / (b*Q), (Q,a) = 1 */
855 : static GEN
856 11829893 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
857 : {
858 11829893 : GEN q = Qdivii(a, b); /* != 0 */
859 11829989 : if (typ(q) == t_INT)
860 : {
861 1538905 : gel(z,1) = gerepileuptoint((pari_sp)Q, q);
862 1538905 : gel(z,2) = Q; return z;
863 : }
864 10291084 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
865 : }
866 : static GEN
867 24655885 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
868 : {
869 24655885 : GEN x1 = gel(x,1), x2 = gel(x,2);
870 24655885 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
871 24655885 : int s = cmpii(x2, y2);
872 :
873 : /* common denominator: (x1 op y1) / x2 */
874 24655886 : if (!s)
875 : {
876 9026700 : pari_sp av = avma;
877 9026700 : return gerepileupto(av, Qdivii(op(x1, y1), x2));
878 : }
879 15629186 : z = cgetg(3, t_FRAC);
880 15629198 : if (s < 0)
881 : {
882 8804523 : Q = dvmdii(y2, x2, &r);
883 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
884 8804427 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
885 2493024 : delta = gcdii(x2,r);
886 : }
887 : else
888 : {
889 6824675 : Q = dvmdii(x2, y2, &r);
890 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
891 6824696 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
892 1306130 : delta = gcdii(y2,r);
893 : }
894 : /* delta = gcd(x2,y2) */
895 3799212 : if (equali1(delta))
896 : { /* numerator is nonzero */
897 1519056 : gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
898 1519057 : gel(z,2) = mulii(x2,y2); return z;
899 : }
900 2280150 : x2 = diviiexact(x2,delta);
901 2280151 : y2 = diviiexact(y2,delta);
902 2280151 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
903 2280152 : q = dvmdii(n, delta, &r);
904 2280150 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
905 2084280 : r = gcdii(delta, r);
906 2084284 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
907 2084284 : return setfrac(z, n, mulii(mulii(x2, y2), delta));
908 : }
909 :
910 : /* assume x2, y2 are t_POLs in the same variable */
911 : static GEN
912 3041742 : add_rfrac(GEN x, GEN y)
913 : {
914 3041742 : pari_sp av = avma;
915 3041742 : GEN x1 = gel(x,1), x2 = gel(x,2);
916 3041742 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
917 :
918 3041742 : delta = RgX_gcd(x2,y2);
919 3041742 : if (!degpol(delta))
920 : {
921 665 : n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
922 665 : d = RgX_mul(x2, y2);
923 665 : return gerepileupto(av, gred_rfrac_simple(n, d));
924 : }
925 3041077 : x2 = RgX_div(x2,delta);
926 3041077 : y2 = RgX_div(y2,delta);
927 3041077 : n = gadd(gmul(x1,y2), gmul(y1,x2));
928 3041077 : if (!signe(n))
929 : {
930 721418 : n = simplify_shallow(n);
931 721418 : if (isexactzero(n))
932 : {
933 721411 : if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
934 7 : return gerepileupto(av, scalarpol(n, varn(x2)));
935 : }
936 7 : return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
937 : }
938 2319659 : if (degpol(n) == 0)
939 1150554 : return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
940 1169105 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
941 1169105 : if (isexactzero(r))
942 : {
943 : GEN z;
944 228010 : d = RgX_mul(x2, y2);
945 : /* "constant" denominator ? */
946 228010 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
947 228010 : return gerepileupto(av, z);
948 : }
949 941095 : r = RgX_gcd(delta, r);
950 941095 : if (degpol(r))
951 : {
952 160726 : n = RgX_div(n, r);
953 160726 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
954 : }
955 : else
956 780369 : d = RgX_mul(gel(x,2), y2);
957 941095 : return gerepileupto(av, gred_rfrac_simple(n, d));
958 : }
959 :
960 : GEN
961 4473397305 : gadd(GEN x, GEN y)
962 : {
963 4473397305 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
964 : pari_sp av, tetpil;
965 : GEN z, p1;
966 :
967 4473397305 : if (tx == ty) switch(tx) /* shortcut to generic case */
968 : {
969 1981321952 : case t_INT: return addii(x,y);
970 1392709627 : case t_REAL: return addrr(x,y);
971 1679202 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
972 1679202 : z = cgetg(3,t_INTMOD);
973 1679202 : if (X==Y || equalii(X,Y))
974 1679188 : return add_intmod_same(z, X, gel(x,2), gel(y,2));
975 14 : gel(z,1) = gcdii(X,Y);
976 14 : warn_coercion(X,Y,gel(z,1));
977 14 : av = avma; p1 = addii(gel(x,2),gel(y,2));
978 14 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
979 : }
980 19667547 : case t_FRAC: return addsub_frac(x,y,addii);
981 270502650 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
982 270511768 : gel(z,2) = gadd(gel(x,2),gel(y,2));
983 270387980 : if (isintzero(gel(z,2)))
984 : {
985 362131 : set_avma((pari_sp)(z+3));
986 362131 : return gadd(gel(x,1),gel(y,1));
987 : }
988 270026315 : gel(z,1) = gadd(gel(x,1),gel(y,1));
989 270024005 : return z;
990 2769673 : case t_PADIC:
991 2769673 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
992 2769463 : return addsub_pp(x,y, addii);
993 672 : case t_QUAD: z = cgetg(4,t_QUAD);
994 672 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
995 672 : gel(z,1) = ZX_copy(gel(x,1));
996 672 : gel(z,2) = gadd(gel(x,2),gel(y,2));
997 672 : gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
998 9857782 : case t_POLMOD:
999 9857782 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1000 9855675 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
1001 2107 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
1002 14006310 : case t_FFELT: return FF_add(x,y);
1003 66605126 : case t_POL:
1004 66605126 : vx = varn(x);
1005 66605126 : vy = varn(y);
1006 66605126 : if (vx != vy) {
1007 872067 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1008 65742 : else return RgX_Rg_add(y, x);
1009 : }
1010 65733059 : return RgX_add(x, y);
1011 3456005 : case t_SER:
1012 3456005 : vx = varn(x);
1013 3456005 : vy = varn(y);
1014 3456005 : if (vx != vy) {
1015 21 : if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1016 21 : else return add_ser_scal(y, x);
1017 : }
1018 3455984 : return ser_add(x, y);
1019 4333629 : case t_RFRAC:
1020 4333629 : vx = varn(gel(x,2));
1021 4333629 : vy = varn(gel(y,2));
1022 4333629 : if (vx != vy) {
1023 1291887 : if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1024 538397 : else return add_rfrac_scal(y, x);
1025 : }
1026 3041742 : return add_rfrac(x,y);
1027 288213 : case t_VEC:
1028 288213 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1029 288213 : return RgV_add(x,y);
1030 1074681 : case t_COL:
1031 1074681 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1032 1074681 : return RgC_add(x,y);
1033 671237 : case t_MAT:
1034 671237 : lx = lg(x);
1035 671237 : if (lg(y) != lx) pari_err_OP("+",x,y);
1036 671237 : if (lx == 1) return cgetg(1, t_MAT);
1037 671237 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1038 671230 : return RgM_add(x,y);
1039 :
1040 0 : default: pari_err_TYPE2("+",x,y);
1041 : }
1042 : /* tx != ty */
1043 706303478 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1044 :
1045 706303478 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1046 : {
1047 626251733 : case t_INT:
1048 : switch(ty)
1049 : {
1050 401842987 : case t_REAL: return addir(x,y);
1051 278105 : case t_INTMOD:
1052 278105 : z = cgetg(3, t_INTMOD);
1053 278105 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1054 43598417 : case t_FRAC: z = cgetg(3,t_FRAC);
1055 43598403 : gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1056 43598393 : gel(z,2) = icopy(gel(y,2)); return z;
1057 174579305 : case t_COMPLEX: return addRc(x, y);
1058 4637450 : case t_PADIC:
1059 4637450 : if (!signe(x)) return gcopy(y);
1060 3487771 : return addQp(x,y);
1061 1113 : case t_QUAD: return addRq(x, y);
1062 1356296 : case t_FFELT: return FF_Z_add(y,x);
1063 : }
1064 :
1065 : case t_REAL:
1066 44896255 : switch(ty)
1067 : {
1068 13118398 : case t_FRAC:
1069 13118398 : if (!signe(gel(y,1))) return rcopy(x);
1070 13118398 : if (!signe(x))
1071 : {
1072 8473 : lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1073 8473 : return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1074 : }
1075 13109925 : av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
1076 13107501 : return gerepile(av,tetpil,divri(z,gel(y,2)));
1077 31777794 : case t_COMPLEX: return addRc(x, y);
1078 63 : case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
1079 :
1080 0 : default: pari_err_TYPE2("+",x,y);
1081 : }
1082 :
1083 17640 : case t_INTMOD:
1084 : switch(ty)
1085 : {
1086 17507 : case t_FRAC: { GEN X = gel(x,1);
1087 17507 : z = cgetg(3, t_INTMOD);
1088 17507 : p1 = Fp_div(gel(y,1), gel(y,2), X);
1089 17507 : return add_intmod_same(z, X, p1, gel(x,2));
1090 : }
1091 14 : case t_FFELT:
1092 14 : if (!equalii(gel(x,1),FF_p_i(y)))
1093 0 : pari_err_OP("+",x,y);
1094 14 : return FF_Z_add(y,gel(x,2));
1095 91 : case t_COMPLEX: return addRc(x, y);
1096 0 : case t_PADIC: { GEN X = gel(x,1);
1097 0 : z = cgetg(3, t_INTMOD);
1098 0 : return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1099 : }
1100 28 : case t_QUAD: return addRq(x, y);
1101 : }
1102 :
1103 : case t_FRAC:
1104 : switch (ty)
1105 : {
1106 10445443 : case t_COMPLEX: return addRc(x, y);
1107 25690 : case t_PADIC:
1108 25690 : if (!signe(gel(x,1))) return gcopy(y);
1109 25690 : return addQp(x,y);
1110 126 : case t_QUAD: return addRq(x, y);
1111 1337 : case t_FFELT: return FF_Q_add(y, x);
1112 : }
1113 :
1114 : case t_FFELT:
1115 0 : pari_err_TYPE2("+",x,y);
1116 :
1117 35 : case t_COMPLEX:
1118 : switch(ty)
1119 : {
1120 28 : case t_PADIC:
1121 28 : return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
1122 7 : case t_QUAD:
1123 7 : lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1124 7 : return gequal0(y)? gcopy(x): addqf(y, x, lx);
1125 : }
1126 :
1127 : case t_PADIC: /* ty == t_QUAD */
1128 0 : return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
1129 : }
1130 : /* tx < ty, !is_const_t(y) */
1131 26388726 : switch(ty)
1132 : {
1133 7993 : case t_MAT:
1134 7993 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1135 7993 : if (isrationalzero(x)) return gcopy(y);
1136 7916 : return RgM_Rg_add(y, x);
1137 202306 : case t_COL:
1138 202306 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1139 202306 : return RgC_Rg_add(y, x);
1140 1838867 : case t_POLMOD: /* is_const_t(tx) in this case */
1141 1838867 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1142 : }
1143 24339560 : if (is_scalar_t(tx)) {
1144 21247240 : if (tx == t_POLMOD)
1145 : {
1146 106400 : vx = varn(gel(x,1));
1147 106400 : vy = gvar(y);
1148 106399 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1149 : else
1150 81304 : if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1151 40124 : return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1152 : }
1153 21140840 : return add_scal(y, x, ty);
1154 : }
1155 : /* x and y are not scalars, ty != t_MAT */
1156 3092290 : vx = gvar(x);
1157 3092290 : vy = gvar(y);
1158 3092290 : if (vx != vy) { /* x or y is treated as a scalar */
1159 22682 : if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1160 32263 : return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1161 32263 : : add_scal(y, x, ty);
1162 : }
1163 : /* vx = vy */
1164 3069608 : switch(tx)
1165 : {
1166 3069146 : case t_POL:
1167 : switch (ty)
1168 : {
1169 5299 : case t_SER:
1170 5299 : if (lg(x) == 2) return gcopy(y);
1171 5285 : i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1172 5285 : i = lg(y) + valp(y) - i;
1173 5285 : if (i < 3) return gcopy(y);
1174 5278 : p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1175 5278 : settyp(p1, t_VECSMALL); /* p1 left on stack */
1176 5278 : return y;
1177 :
1178 3063847 : case t_RFRAC: return add_rfrac_scal(y, x);
1179 : }
1180 0 : break;
1181 :
1182 462 : case t_SER:
1183 462 : if (ty == t_RFRAC)
1184 : {
1185 : long vn, vd;
1186 462 : av = avma;
1187 462 : vn = gval(gel(y,1), vy);
1188 462 : vd = RgX_valrem_inexact(gel(y,2), NULL);
1189 462 : if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
1190 :
1191 462 : l = lg(x) + valp(x) - (vn - vd);
1192 462 : if (l < 3) { set_avma(av); return gcopy(x); }
1193 462 : return gerepileupto(av, gadd(x, rfrac_to_ser_i(y, l)));
1194 : }
1195 0 : break;
1196 : }
1197 0 : pari_err_TYPE2("+",x,y);
1198 : return NULL; /* LCOV_EXCL_LINE */
1199 : }
1200 :
1201 : GEN
1202 283926422 : gaddsg(long x, GEN y)
1203 : {
1204 283926422 : long ty = typ(y);
1205 : GEN z;
1206 :
1207 283926422 : switch(ty)
1208 : {
1209 133240627 : case t_INT: return addsi(x,y);
1210 126998568 : case t_REAL: return addsr(x,y);
1211 14 : case t_INTMOD:
1212 14 : z = cgetg(3, t_INTMOD);
1213 14 : return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1214 14710451 : case t_FRAC: z = cgetg(3,t_FRAC);
1215 14710451 : gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1216 14710451 : gel(z,2) = icopy(gel(y,2)); return z;
1217 8235120 : case t_COMPLEX:
1218 8235120 : z = cgetg(3, t_COMPLEX);
1219 8235120 : gel(z,1) = gaddsg(x, gel(y,1));
1220 8235120 : gel(z,2) = gcopy(gel(y,2)); return z;
1221 :
1222 741642 : default: return gadd(stoi(x), y);
1223 : }
1224 : }
1225 :
1226 : GEN
1227 3144202 : gsubsg(long x, GEN y)
1228 : {
1229 : GEN z, a, b;
1230 : pari_sp av;
1231 :
1232 3144202 : switch(typ(y))
1233 : {
1234 274342 : case t_INT: return subsi(x,y);
1235 1212608 : case t_REAL: return subsr(x,y);
1236 56 : case t_INTMOD:
1237 56 : z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1238 56 : return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1239 729137 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1240 729137 : gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
1241 729134 : gel(z,2) = icopy(gel(y,2)); return z;
1242 891337 : case t_COMPLEX:
1243 891337 : z = cgetg(3, t_COMPLEX);
1244 891337 : gel(z,1) = gsubsg(x, gel(y,1));
1245 891337 : gel(z,2) = gneg(gel(y,2)); return z;
1246 : }
1247 36722 : av = avma;
1248 36722 : return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
1249 : }
1250 :
1251 : /********************************************************************/
1252 : /** **/
1253 : /** SUBTRACTION **/
1254 : /** **/
1255 : /********************************************************************/
1256 :
1257 : GEN
1258 3427362920 : gsub(GEN x, GEN y)
1259 : {
1260 3427362920 : long tx = typ(x), ty = typ(y);
1261 : pari_sp av;
1262 : GEN z;
1263 3427362920 : if (tx == ty) switch(tx) /* shortcut to generic case */
1264 : {
1265 2743228627 : case t_INT: return subii(x,y);
1266 495630649 : case t_REAL: return subrr(x,y);
1267 1153268 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1268 1153268 : z = cgetg(3,t_INTMOD);
1269 1153268 : if (X==Y || equalii(X,Y))
1270 1153254 : return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1271 14 : gel(z,1) = gcdii(X,Y);
1272 14 : warn_coercion(X,Y,gel(z,1));
1273 14 : av = avma; p1 = subii(gel(x,2),gel(y,2));
1274 14 : gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
1275 : }
1276 4988368 : case t_FRAC: return addsub_frac(x,y, subii);
1277 107628951 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1278 107669373 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1279 107558093 : if (isintzero(gel(z,2)))
1280 : {
1281 21237 : set_avma((pari_sp)(z+3));
1282 21237 : return gsub(gel(x,1),gel(y,1));
1283 : }
1284 107540652 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1285 107526884 : return z;
1286 7762227 : case t_PADIC:
1287 7762227 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
1288 7762227 : return addsub_pp(x,y, subii);
1289 168 : case t_QUAD: z = cgetg(4,t_QUAD);
1290 168 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1291 168 : gel(z,1) = ZX_copy(gel(x,1));
1292 168 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1293 168 : gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1294 793265 : case t_POLMOD:
1295 793265 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1296 793230 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1297 35 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1298 203128 : case t_FFELT: return FF_sub(x,y);
1299 10219580 : case t_POL: {
1300 10219580 : long vx = varn(x);
1301 10219580 : long vy = varn(y);
1302 10219580 : if (vx != vy) {
1303 30219 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1304 5614 : else return Rg_RgX_sub(x, y);
1305 : }
1306 10189361 : return RgX_sub(x, y);
1307 : }
1308 297384 : case t_VEC:
1309 297384 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1310 297384 : return RgV_sub(x,y);
1311 6534220 : case t_COL:
1312 6534220 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1313 6534220 : return RgC_sub(x,y);
1314 234822 : case t_MAT: {
1315 234822 : long lx = lg(x);
1316 234822 : if (lg(y) != lx) pari_err_OP("+",x,y);
1317 234821 : if (lx == 1) return cgetg(1, t_MAT);
1318 155735 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1319 155735 : return RgM_sub(x,y);
1320 : }
1321 2452440 : case t_RFRAC: case t_SER: break;
1322 :
1323 0 : default: pari_err_TYPE2("+",x,y);
1324 : }
1325 47891231 : av = avma;
1326 50343671 : return gerepileupto(av, gadd(x,gneg_i(y)));
1327 : }
1328 :
1329 : /********************************************************************/
1330 : /** **/
1331 : /** MULTIPLICATION **/
1332 : /** **/
1333 : /********************************************************************/
1334 : static GEN
1335 302933 : mul_ser_scal(GEN y, GEN x)
1336 : {
1337 : long l, i;
1338 : GEN z;
1339 302933 : if (isexactzero(x)) return gmul(Rg_get_0(y), x);
1340 299685 : if (ser_isexactzero(y))
1341 : {
1342 476 : z = scalarser(lg(y) == 2? Rg_get_0(x): gmul(gel(y,2), x), varn(y), 1);
1343 476 : setvalp(z, valp(y)); return z;
1344 : }
1345 299209 : z = cgetg_copy(y, &l); z[1] = y[1];
1346 4167984 : for (i = 2; i < l; i++) gel(z,i) = gmul(gel(y,i), x);
1347 299209 : return normalizeser(z);
1348 : }
1349 : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1350 : * [n/d a valid RFRAC] */
1351 : static GEN
1352 10451122 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1353 : {
1354 10451122 : pari_sp av = avma;
1355 : GEN z;
1356 :
1357 10451122 : switch(typ(x))
1358 : {
1359 21 : case t_PADIC:
1360 21 : n = gmul(n, x);
1361 21 : d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
1362 21 : return gerepileupto(av, gdiv(n,d));
1363 :
1364 490 : case t_INTMOD: case t_POLMOD:
1365 490 : n = gmul(n, x);
1366 490 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1367 490 : return gerepileupto(av, gdiv(n,d));
1368 : }
1369 10450611 : z = gred_rfrac2(x, d);
1370 10450611 : n = simplify_shallow(n);
1371 10450611 : if (typ(z) == t_RFRAC)
1372 : {
1373 7933879 : n = gmul(gel(z,1), n);
1374 7933879 : d = gel(z,2);
1375 7933879 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1376 0 : z = RgX_Rg_div(n, d);
1377 : else
1378 7933879 : z = gred_rfrac_simple(n, d);
1379 : }
1380 : else
1381 2516732 : z = gmul(z, n);
1382 10450611 : return gerepileupto(av, z);
1383 : }
1384 : static GEN
1385 111560892 : mul_scal(GEN y, GEN x, long ty)
1386 : {
1387 111560892 : switch(ty)
1388 : {
1389 102679427 : case t_POL:
1390 102679427 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1391 99527959 : return RgX_Rg_mul(y, x);
1392 294848 : case t_SER: return mul_ser_scal(y, x);
1393 8586627 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1394 14 : case t_QFB:
1395 14 : if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1396 : }
1397 12 : pari_err_TYPE2("*",x,y);
1398 : return NULL; /* LCOV_EXCL_LINE */
1399 : }
1400 :
1401 : static GEN
1402 160865 : mul_gen_rfrac(GEN X, GEN Y)
1403 : {
1404 160865 : GEN y1 = gel(Y,1), y2 = gel(Y,2);
1405 160865 : long vx = gvar(X), vy = varn(y2);
1406 166626 : return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1407 5761 : gred_rfrac_simple(gmul(y1,X), y2);
1408 : }
1409 : /* (x1/x2) * (y1/y2) */
1410 : static GEN
1411 7907093 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1412 : {
1413 : GEN z, X, Y;
1414 7907093 : pari_sp av = avma;
1415 :
1416 7907093 : X = gred_rfrac2(x1, y2);
1417 7907093 : Y = gred_rfrac2(y1, x2);
1418 7907093 : if (typ(X) == t_RFRAC)
1419 : {
1420 6628032 : if (typ(Y) == t_RFRAC) {
1421 6562456 : x1 = gel(X,1);
1422 6562456 : x2 = gel(X,2);
1423 6562456 : y1 = gel(Y,1);
1424 6562456 : y2 = gel(Y,2);
1425 6562456 : z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1426 : } else
1427 65576 : z = mul_gen_rfrac(Y, X);
1428 : }
1429 1279061 : else if (typ(Y) == t_RFRAC)
1430 95289 : z = mul_gen_rfrac(X, Y);
1431 : else
1432 1183772 : z = gmul(X, Y);
1433 7907093 : return gerepileupto(av, z);
1434 : }
1435 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1436 : static GEN
1437 269687 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1438 : {
1439 269687 : pari_sp av = avma;
1440 269687 : GEN X = gred_rfrac2(x1, y2);
1441 269687 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1442 : {
1443 263121 : x2 = RgX_mul(gel(X,2), x2);
1444 263121 : x1 = gel(X,1);
1445 : }
1446 : else
1447 6566 : x1 = X;
1448 269687 : return gerepileupto(av, gred_rfrac_simple(x1, x2));
1449 : }
1450 :
1451 : /* Mod(y, Y) * x, assuming x scalar */
1452 : static GEN
1453 2571802 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1454 2571802 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1455 :
1456 : /* cf mulqq */
1457 : static GEN
1458 5854148 : quad_polmod_mul(GEN P, GEN x, GEN y)
1459 : {
1460 5854148 : GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
1461 5854148 : pari_sp tetpil, av = avma;
1462 5854148 : T[1] = x[1];
1463 5854148 : p2 = gmul(gel(x,2), gel(y,2));
1464 5854148 : p3 = gmul(gel(x,3), gel(y,3));
1465 5854148 : p1 = gmul(gneg_i(c),p3);
1466 : /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
1467 5854148 : if (typ(b) == t_INT)
1468 : {
1469 5854127 : if (signe(b))
1470 : {
1471 4284657 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1472 4284657 : if (is_pm1(b))
1473 : {
1474 4283999 : if (signe(b) > 0) p3 = gneg(p3);
1475 : }
1476 : else
1477 658 : p3 = gmul(negi(b), p3);
1478 : }
1479 : else
1480 : {
1481 1569470 : p3 = gmul(gel(x,2),gel(y,3));
1482 1569470 : p4 = gmul(gel(x,3),gel(y,2));
1483 : }
1484 : }
1485 : else
1486 : {
1487 21 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1488 21 : p3 = gmul(gneg_i(b), p3);
1489 : }
1490 5854148 : tetpil = avma;
1491 5854148 : gel(T,2) = gadd(p2, p1);
1492 5854148 : gel(T,3) = gadd(p4, p3);
1493 5854148 : gerepilecoeffssp(av,tetpil,T+2,2);
1494 5854148 : return normalizepol_lg(T,4);
1495 : }
1496 : /* Mod(x,T) * Mod(y,T) */
1497 : static GEN
1498 7782603 : mul_polmod_same(GEN T, GEN x, GEN y)
1499 : {
1500 7782603 : GEN z = cgetg(3,t_POLMOD), a;
1501 7782603 : long v = varn(T), lx = lg(x), ly = lg(y);
1502 7782603 : gel(z,1) = RgX_copy(T);
1503 : /* x * y mod T optimised */
1504 7782603 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1505 7163285 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1506 1489363 : a = gmul(x, y);
1507 : else
1508 : {
1509 6293240 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1510 5849003 : a = quad_polmod_mul(T, x, y);
1511 : else
1512 444237 : a = RgXQ_mul(x, y, gel(z,1));
1513 : }
1514 7782603 : gel(z,2) = a; return z;
1515 : }
1516 : static GEN
1517 46880 : sqr_polmod(GEN T, GEN x)
1518 : {
1519 46880 : GEN a, z = cgetg(3,t_POLMOD);
1520 46880 : gel(z,1) = RgX_copy(T);
1521 46880 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1522 14226 : a = gsqr(x);
1523 : else
1524 : {
1525 32654 : pari_sp av = avma;
1526 32654 : a = RgXQ_sqr(x, gel(z,1));
1527 32654 : a = gerepileupto(av, a);
1528 : }
1529 46880 : gel(z,2) = a; return z;
1530 : }
1531 : /* Mod(x,X) * Mod(y,Y) */
1532 : static GEN
1533 2678207 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1534 : {
1535 2678207 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1536 2678207 : long vx = varn(X), vy = varn(Y);
1537 2678207 : GEN z = cgetg(3,t_POLMOD);
1538 :
1539 2678207 : if (vx==vy) {
1540 : pari_sp av;
1541 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
1542 14 : warn_coercion(X,Y,gel(z,1));
1543 14 : gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
1544 14 : return z;
1545 : }
1546 2678193 : if (varncmp(vx, vy) < 0)
1547 414316 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1548 : else
1549 2263877 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1550 2678193 : gel(z,2) = gmul(x, y); return z;
1551 : }
1552 :
1553 : #if 0 /* used by 3M only */
1554 : /* set z = x+y and return 1 if x,y have the same sign
1555 : * set z = x-y and return 0 otherwise */
1556 : static int
1557 : did_add(GEN x, GEN y, GEN *z)
1558 : {
1559 : long tx = typ(x), ty = typ(y);
1560 : if (tx == ty) switch(tx)
1561 : {
1562 : case t_INT: *z = addii(x,y); return 1;
1563 : case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1564 : case t_REAL:
1565 : if (signe(x) == -signe(y))
1566 : { *z = subrr(x,y); return 0; }
1567 : else
1568 : { *z = addrr(x,y); return 1; }
1569 : }
1570 : if (tx == t_REAL) switch(ty)
1571 : {
1572 : case t_INT:
1573 : if (signe(x) == -signe(y))
1574 : { *z = subri(x,y); return 0; }
1575 : else
1576 : { *z = addri(x,y); return 1; }
1577 : case t_FRAC:
1578 : if (signe(x) == -signe(gel(y,1)))
1579 : { *z = gsub(x,y); return 0; }
1580 : else
1581 : { *z = gadd(x,y); return 1; }
1582 : }
1583 : else if (ty == t_REAL) switch(tx)
1584 : {
1585 : case t_INT:
1586 : if (signe(x) == -signe(y))
1587 : { *z = subir(x,y); return 0; }
1588 : else
1589 : { *z = addir(x,y); return 1; }
1590 : case t_FRAC:
1591 : if (signe(gel(x,1)) == -signe(y))
1592 : { *z = gsub(x,y); return 0; }
1593 : else
1594 : { *z = gadd(x,y); return 1; }
1595 : }
1596 : *z = gadd(x,y); return 1;
1597 : }
1598 : #endif
1599 : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1600 : static GEN
1601 10995999 : mulcIR(GEN x, GEN y)
1602 : {
1603 10995999 : GEN z = cgetg(3,t_COMPLEX);
1604 10996196 : pari_sp av = avma;
1605 10996196 : gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
1606 10996414 : gel(z,2) = gmul(y, gel(x,1));
1607 10996029 : return z;
1608 :
1609 : }
1610 : /* x,y COMPLEX */
1611 : static GEN
1612 246951223 : mulcc(GEN x, GEN y)
1613 : {
1614 246951223 : GEN xr = gel(x,1), xi = gel(x,2);
1615 246951223 : GEN yr = gel(y,1), yi = gel(y,2);
1616 : GEN p1, p2, p3, p4, z;
1617 : pari_sp tetpil, av;
1618 :
1619 246951223 : if (isintzero(xr))
1620 : {
1621 15315067 : if (isintzero(yr)) {
1622 7160003 : av = avma;
1623 7160003 : return gerepileupto(av, gneg(gmul(xi,yi)));
1624 : }
1625 8154855 : return mulcIR(y, xi);
1626 : }
1627 231636687 : if (isintzero(yr)) return mulcIR(x, yi);
1628 :
1629 228805263 : z = cgetg(3,t_COMPLEX); av = avma;
1630 : #if 0
1631 : /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1632 : * e.g. xr + xi if exponents differ */
1633 : if (did_add(xr, xi, &p3))
1634 : {
1635 : if (did_add(yr, yi, &p4)) {
1636 : /* R = xr*yr - xi*yi
1637 : * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1638 : p1 = gmul(xr,yr);
1639 : p2 = gmul(xi,yi); p2 = gneg(p2);
1640 : p3 = gmul(p3, p4);
1641 : p4 = gsub(p2, p1);
1642 : } else {
1643 : /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1644 : * I = xr*yi + xi*yr */
1645 : p1 = gmul(p3,p4);
1646 : p3 = gmul(xr,yi);
1647 : p4 = gmul(xi,yr);
1648 : p2 = gsub(p3, p4);
1649 : }
1650 : } else {
1651 : if (did_add(yr, yi, &p4)) {
1652 : /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1653 : * I = xr*yi +xi*yr */
1654 : p1 = gmul(p3,p4);
1655 : p3 = gmul(xr,yi);
1656 : p4 = gmul(xi,yr);
1657 : p2 = gsub(p4, p3);
1658 : } else {
1659 : /* R = xr*yr - xi*yi
1660 : * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1661 : p3 = gneg( gmul(p3, p4) );
1662 : p1 = gmul(xr,yr);
1663 : p2 = gmul(xi,yi);
1664 : p4 = gadd(p1, p2);
1665 :
1666 : p2 = gneg(p2);
1667 : }
1668 : }
1669 : tetpil = avma;
1670 : gel(z,1) = gadd(p1,p2);
1671 : gel(z,2) = gadd(p3,p4);
1672 : #else
1673 228826886 : if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1674 : { /* 3M formula */
1675 557442 : p3 = addii(xr,xi);
1676 557442 : p4 = addii(yr,yi);
1677 557442 : p1 = mulii(xr,yr);
1678 557442 : p2 = mulii(xi,yi);
1679 557442 : p3 = mulii(p3,p4);
1680 557442 : p4 = addii(p2,p1);
1681 557442 : tetpil = avma;
1682 557442 : gel(z,1) = subii(p1,p2);
1683 557442 : gel(z,2) = subii(p3,p4);
1684 557442 : if (!signe(gel(z,2)))
1685 113057 : return gerepileuptoint((pari_sp)(z+3), gel(z,1));
1686 : }
1687 : else
1688 : { /* naive 4M formula: avoid all loss of accuracy */
1689 228269444 : p1 = gmul(xr,yr);
1690 228196583 : p2 = gmul(xi,yi);
1691 228176803 : p3 = gmul(xr,yi);
1692 228180222 : p4 = gmul(xi,yr);
1693 228188600 : tetpil = avma;
1694 228188600 : gel(z,1) = gsub(p1,p2);
1695 228098981 : gel(z,2) = gadd(p3,p4);
1696 228104692 : if (isintzero(gel(z,2)))
1697 : {
1698 31387 : cgiv(gel(z,2));
1699 31387 : return gerepileupto((pari_sp)(z+3), gel(z,1));
1700 : }
1701 : }
1702 : #endif
1703 :
1704 228519035 : gerepilecoeffssp(av,tetpil, z+1,2); return z;
1705 : }
1706 : /* x,y PADIC */
1707 : static GEN
1708 16846772 : mulpp(GEN x, GEN y) {
1709 16846772 : long l = valp(x) + valp(y);
1710 : pari_sp av;
1711 : GEN z, t;
1712 16846772 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
1713 16846730 : if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
1714 16630093 : if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
1715 :
1716 16293128 : t = (precp(x) > precp(y))? y: x;
1717 16293128 : z = cgetp(t); setvalp(z,l); av = avma;
1718 16293231 : affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
1719 16292991 : set_avma(av); return z;
1720 : }
1721 : /* x,y QUAD */
1722 : static GEN
1723 1106 : mulqq(GEN x, GEN y) {
1724 1106 : GEN z = cgetg(4,t_QUAD);
1725 1106 : GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
1726 : pari_sp av, tetpil;
1727 1106 : if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
1728 :
1729 1106 : gel(z,1) = ZX_copy(P); av = avma;
1730 1106 : p2 = gmul(gel(x,2),gel(y,2));
1731 1106 : p3 = gmul(gel(x,3),gel(y,3));
1732 1106 : p1 = gmul(gneg_i(c),p3);
1733 :
1734 1106 : if (signe(b))
1735 987 : p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
1736 : else
1737 : {
1738 119 : p3 = gmul(gel(x,2),gel(y,3));
1739 119 : p4 = gmul(gel(x,3),gel(y,2));
1740 : }
1741 1106 : tetpil = avma;
1742 1106 : gel(z,2) = gadd(p2,p1);
1743 1106 : gel(z,3) = gadd(p4,p3);
1744 1106 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
1745 : }
1746 :
1747 : GEN
1748 10320816 : mulcxI(GEN x)
1749 : {
1750 : GEN z;
1751 10320816 : switch(typ(x))
1752 : {
1753 1561452 : case t_INT: case t_REAL: case t_FRAC:
1754 1561452 : return mkcomplex(gen_0, x);
1755 8759908 : case t_COMPLEX:
1756 8759908 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1757 8759278 : z = cgetg(3,t_COMPLEX);
1758 8759846 : gel(z,1) = gneg(gel(x,2));
1759 8760331 : gel(z,2) = gel(x,1); return z;
1760 48 : default:
1761 48 : return gmul(gen_I(), x);
1762 : }
1763 : }
1764 : GEN
1765 196014 : mulcxmI(GEN x)
1766 : {
1767 : GEN z;
1768 196014 : switch(typ(x))
1769 : {
1770 84922 : case t_INT: case t_REAL: case t_FRAC:
1771 84922 : return mkcomplex(gen_0, gneg(x));
1772 111029 : case t_COMPLEX:
1773 111029 : if (isintzero(gel(x,1))) return gel(x,2);
1774 108418 : z = cgetg(3,t_COMPLEX);
1775 108418 : gel(z,1) = gel(x,2);
1776 108418 : gel(z,2) = gneg(gel(x,1)); return z;
1777 63 : default:
1778 63 : return gmul(mkcomplex(gen_0, gen_m1), x);
1779 : }
1780 : }
1781 : /* x * I^k */
1782 : GEN
1783 259031 : mulcxpowIs(GEN x, long k)
1784 : {
1785 259031 : switch (k & 3)
1786 : {
1787 49649 : case 1: return mulcxI(x);
1788 46337 : case 2: return gneg(x);
1789 49968 : case 3: return mulcxmI(x);
1790 : }
1791 113077 : return x;
1792 : }
1793 :
1794 : static GEN
1795 5394735 : init_ser(long l, long v, long e)
1796 : {
1797 5394735 : GEN z = cgetg(l, t_SER);
1798 5394735 : z[1] = evalvalp(e) | evalvarn(v) | evalsigne(1); return z;
1799 : }
1800 :
1801 : /* fill in coefficients of t_SER z from coeffs of t_POL y */
1802 : static GEN
1803 5379077 : fill_ser(GEN z, GEN y)
1804 : {
1805 5379077 : long i, lx = lg(z), ly = lg(y);
1806 5379077 : if (ly >= lx) {
1807 24020714 : for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
1808 : } else {
1809 328331 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1810 237900 : for ( ; i < lx; i++) gel(z,i) = gen_0;
1811 : }
1812 5379077 : return normalizeser(z);
1813 : }
1814 :
1815 : GEN
1816 7748018530 : gmul(GEN x, GEN y)
1817 : {
1818 : long tx, ty, lx, ly, vx, vy, i, l;
1819 : pari_sp av, tetpil;
1820 : GEN z, p1;
1821 :
1822 7748018530 : if (x == y) return gsqr(x);
1823 6875357121 : tx = typ(x); ty = typ(y);
1824 6875357121 : if (tx == ty) switch(tx)
1825 : {
1826 3418375397 : case t_INT: return mulii(x,y);
1827 1632440238 : case t_REAL: return mulrr(x,y);
1828 2048342 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1829 2048342 : z = cgetg(3,t_INTMOD);
1830 2048342 : if (X==Y || equalii(X,Y))
1831 2048306 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1832 35 : gel(z,1) = gcdii(X,Y);
1833 35 : warn_coercion(X,Y,gel(z,1));
1834 35 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1835 35 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
1836 : }
1837 20622238 : case t_FRAC:
1838 : {
1839 20622238 : GEN x1 = gel(x,1), x2 = gel(x,2);
1840 20622238 : GEN y1 = gel(y,1), y2 = gel(y,2);
1841 20622238 : z=cgetg(3,t_FRAC);
1842 20622238 : p1 = gcdii(x1, y2);
1843 20622229 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1844 20622228 : p1 = gcdii(x2, y1);
1845 20622230 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1846 20622229 : tetpil = avma;
1847 20622229 : gel(z,2) = mulii(x2,y2);
1848 20622225 : gel(z,1) = mulii(x1,y1);
1849 20622225 : fix_frac_if_int_GC(z,tetpil); return z;
1850 : }
1851 239363270 : case t_COMPLEX: return mulcc(x, y);
1852 10344835 : case t_PADIC: return mulpp(x, y);
1853 882 : case t_QUAD: return mulqq(x, y);
1854 14523401 : case t_FFELT: return FF_mul(x, y);
1855 10213227 : case t_POLMOD:
1856 10213227 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1857 7535020 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1858 2678207 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1859 43948738 : case t_POL:
1860 43948738 : vx = varn(x);
1861 43948738 : vy = varn(y);
1862 43948738 : if (vx != vy) {
1863 4867710 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1864 2082600 : else return RgX_Rg_mul(y, x);
1865 : }
1866 39081028 : return RgX_mul(x, y);
1867 :
1868 4224360 : case t_SER: {
1869 4224360 : vx = varn(x);
1870 4224360 : vy = varn(y);
1871 4224360 : if (vx != vy) {
1872 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1873 3675 : return mul_ser_scal(y, x);
1874 : }
1875 4220685 : lx = minss(lg(x), lg(y));
1876 4220685 : if (lx == 2) return zeroser(vx, valp(x)+valp(y));
1877 4220671 : av = avma; z = init_ser(lx, vx, valp(x)+valp(y));
1878 4220671 : x = ser2pol_i(x, lx);
1879 4220671 : y = ser2pol_i(y, lx);
1880 4220671 : y = RgXn_mul(x, y, lx-2);
1881 4220671 : return gerepilecopy(av, fill_ser(z,y));
1882 : }
1883 55815 : case t_VEC: /* handle extended t_QFB */
1884 55815 : case t_QFB: return qfbcomp(x,y);
1885 6720458 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1886 2749259 : case t_MAT: return RgM_mul(x, y);
1887 :
1888 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1889 1631 : z = cgetg_copy(x, &l);
1890 1631 : if (l != lg(y)) break;
1891 17325 : for (i=1; i<l; i++)
1892 : {
1893 15694 : long yi = y[i];
1894 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1895 15694 : z[i] = x[yi];
1896 : }
1897 1631 : return z;
1898 :
1899 0 : default:
1900 0 : pari_err_TYPE2("*",x,y);
1901 : }
1902 : /* tx != ty */
1903 1471528078 : if (is_const_t(ty) && is_const_t(tx)) {
1904 1349497520 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1905 1349497520 : switch(tx) {
1906 1243608215 : case t_INT:
1907 : switch(ty)
1908 : {
1909 873849758 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1910 1169776 : case t_INTMOD:
1911 1169776 : z = cgetg(3, t_INTMOD);
1912 1169775 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1913 62194224 : case t_FRAC:
1914 62194224 : if (!signe(x)) return gen_0;
1915 44976705 : z=cgetg(3,t_FRAC);
1916 44976891 : p1 = gcdii(x,gel(y,2));
1917 44975585 : if (equali1(p1))
1918 : {
1919 27376779 : set_avma((pari_sp)z);
1920 27376809 : gel(z,2) = icopy(gel(y,2));
1921 27376951 : gel(z,1) = mulii(gel(y,1), x);
1922 : }
1923 : else
1924 : {
1925 17599156 : x = diviiexact(x,p1); tetpil = avma;
1926 17598424 : gel(z,2) = diviiexact(gel(y,2), p1);
1927 17598294 : gel(z,1) = mulii(gel(y,1), x);
1928 17598714 : fix_frac_if_int_GC(z,tetpil);
1929 : }
1930 44976416 : return z;
1931 299284639 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1932 5671998 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1933 1701 : case t_QUAD: return mulRq(x,y);
1934 1622167 : case t_FFELT: return FF_Z_mul(y,x);
1935 : }
1936 :
1937 : case t_REAL:
1938 97616005 : switch(ty)
1939 : {
1940 29262858 : case t_FRAC: return mulrfrac(x, y);
1941 68353105 : case t_COMPLEX: return mulRc(x, y);
1942 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1943 21 : default: pari_err_TYPE2("*",x,y);
1944 : }
1945 :
1946 8064 : case t_INTMOD:
1947 : switch(ty)
1948 : {
1949 7350 : case t_FRAC: { GEN X = gel(x,1);
1950 7350 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1951 7350 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1952 : }
1953 49 : case t_COMPLEX: return mulRc_direct(x,y);
1954 252 : case t_PADIC: { GEN X = gel(x,1);
1955 252 : z = cgetg(3, t_INTMOD);
1956 252 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1957 : }
1958 63 : case t_QUAD: return mulRq(x, y);
1959 350 : case t_FFELT:
1960 350 : if (!equalii(gel(x,1),FF_p_i(y)))
1961 0 : pari_err_OP("*",x,y);
1962 350 : return FF_Z_mul(y,gel(x,2));
1963 : }
1964 :
1965 : case t_FRAC:
1966 : switch(ty)
1967 : {
1968 5103845 : case t_COMPLEX: return mulRc(x, y);
1969 3335979 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
1970 343 : case t_QUAD: return mulRq(x, y);
1971 1995 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
1972 : }
1973 :
1974 : case t_FFELT:
1975 35 : pari_err_TYPE2("*",x,y);
1976 :
1977 21 : case t_COMPLEX:
1978 : switch(ty)
1979 : {
1980 14 : case t_PADIC:
1981 14 : return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
1982 7 : case t_QUAD:
1983 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
1984 7 : return mulqf(y, x, lx);
1985 : }
1986 :
1987 : case t_PADIC: /* ty == t_QUAD */
1988 28 : return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
1989 : }
1990 : }
1991 :
1992 125740635 : if (is_matvec_t(ty))
1993 : {
1994 9001154 : if (!is_matvec_t(tx))
1995 : {
1996 8658933 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
1997 8658931 : z = cgetg_copy(y, &ly);
1998 1407306797 : for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
1999 8658471 : return z;
2000 : }
2001 342220 : switch(tx)
2002 : {
2003 120754 : case t_VEC:
2004 : switch(ty) {
2005 15456 : case t_COL: return RgV_RgC_mul(x,y);
2006 105298 : case t_MAT: return RgV_RgM_mul(x,y);
2007 : }
2008 0 : break;
2009 1687 : case t_COL:
2010 : switch(ty) {
2011 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2012 0 : case t_MAT: return RgC_RgM_mul(x,y);
2013 : }
2014 0 : break;
2015 219796 : case t_MAT:
2016 : switch(ty) {
2017 0 : case t_VEC: return RgM_RgV_mul(x,y);
2018 219796 : case t_COL: return RgM_RgC_mul(x,y);
2019 : }
2020 : }
2021 : }
2022 117205253 : if (is_matvec_t(tx))
2023 : {
2024 1109328 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2025 1109542 : z = cgetg_copy(x, &lx);
2026 5290268 : for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
2027 1109249 : return z;
2028 : }
2029 116095979 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2030 : /* tx < ty, !ismatvec(x and y) */
2031 :
2032 116095979 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2033 2524433 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2034 113571546 : if (is_scalar_t(tx)) {
2035 110281440 : if (tx == t_POLMOD) {
2036 3098711 : vx = varn(gel(x,1));
2037 3098711 : vy = gvar(y);
2038 3098711 : if (vx != vy) {
2039 2851688 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2040 47369 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2041 : }
2042 : /* error if ty == t_SER */
2043 247023 : av = avma; y = gmod(y, gel(x,1));
2044 247016 : return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2045 : }
2046 107182729 : return mul_scal(y, x, ty);
2047 : }
2048 :
2049 : /* x and y are not scalars, nor matvec */
2050 3289881 : vx = gvar(x);
2051 3289914 : vy = gvar(y);
2052 3289914 : if (vx != vy) /* x or y is treated as a scalar */
2053 2780047 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2054 2780047 : : mul_scal(y, x, ty);
2055 : /* vx = vy */
2056 1871138 : switch(tx)
2057 : {
2058 1871110 : case t_POL:
2059 : switch (ty)
2060 : {
2061 6685 : case t_SER:
2062 : {
2063 : long v;
2064 6685 : av = avma; v = RgX_valrem(x, &x);
2065 6685 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2066 6671 : v += valp(y); ly = lg(y);
2067 6671 : if (ly == 2) { set_avma(av); return zeroser(vx, v); }
2068 6440 : if (degpol(x))
2069 : {
2070 2030 : z = init_ser(ly, vx, v);
2071 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2072 2030 : return gerepilecopy(av, fill_ser(z, x));
2073 : }
2074 : /* take advantage of x = c*t^v */
2075 4410 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2076 4410 : setvalp(y, v); return y;
2077 : }
2078 :
2079 1864425 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2080 : }
2081 0 : break;
2082 :
2083 28 : case t_SER:
2084 : switch (ty)
2085 : {
2086 28 : case t_RFRAC:
2087 28 : av = avma;
2088 28 : return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2089 : }
2090 0 : break;
2091 : }
2092 0 : pari_err_TYPE2("*",x,y);
2093 : return NULL; /* LCOV_EXCL_LINE */
2094 : }
2095 :
2096 : /* return a nonnormalized result */
2097 : GEN
2098 112829 : sqr_ser_part(GEN x, long l1, long l2)
2099 : {
2100 : long i, j, l;
2101 : pari_sp av;
2102 : GEN Z, z, p1, p2;
2103 : long mi;
2104 112829 : if (l2 < l1) return zeroser(varn(x), 2*valp(x));
2105 112815 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2106 112815 : Z = cgetg(l2-l1+3, t_SER);
2107 112815 : Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
2108 112815 : z = Z + 2-l1;
2109 112815 : x += 2; mi = 0;
2110 441896 : for (i=0; i<l1; i++)
2111 : {
2112 329081 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2113 : }
2114 :
2115 712344 : for (i=l1; i<=l2; i++)
2116 : {
2117 599529 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2118 599529 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2119 1248434 : for (j=i-mi; j<=minss(l,mi); j++)
2120 648905 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2121 599529 : p1 = gshift(p1,1);
2122 599529 : if ((i&1) == 0 && p2[i>>1])
2123 88356 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2124 599529 : gel(z,i) = gerepileupto(av,p1);
2125 : }
2126 112815 : return Z;
2127 : }
2128 :
2129 : GEN
2130 1055889643 : gsqr(GEN x)
2131 : {
2132 : long i, lx;
2133 : pari_sp av, tetpil;
2134 : GEN z, p1, p2, p3, p4;
2135 :
2136 1055889643 : switch(typ(x))
2137 : {
2138 882119462 : case t_INT: return sqri(x);
2139 160220509 : case t_REAL: return sqrr(x);
2140 142692 : case t_INTMOD: { GEN X = gel(x,1);
2141 142692 : z = cgetg(3,t_INTMOD);
2142 142693 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
2143 142694 : gel(z,1) = icopy(X); return z;
2144 : }
2145 2976002 : case t_FRAC: return sqrfrac(x);
2146 :
2147 6592751 : case t_COMPLEX:
2148 6592751 : if (isintzero(gel(x,1))) {
2149 209263 : av = avma;
2150 209263 : return gerepileupto(av, gneg(gsqr(gel(x,2))));
2151 : }
2152 6383493 : z = cgetg(3,t_COMPLEX); av = avma;
2153 6383482 : p1 = gadd(gel(x,1),gel(x,2));
2154 6383321 : p2 = gsub(gel(x,1), gel(x,2));
2155 6383245 : p3 = gmul(gel(x,1),gel(x,2));
2156 6383347 : tetpil = avma;
2157 6383347 : gel(z,1) = gmul(p1,p2);
2158 6383370 : gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
2159 :
2160 2583 : case t_PADIC:
2161 2583 : z = cgetg(5,t_PADIC);
2162 2583 : i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
2163 2583 : if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
2164 2583 : z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
2165 2583 : gel(z,2) = icopy(gel(x,2));
2166 2583 : gel(z,3) = shifti(gel(x,3), i); av = avma;
2167 2583 : gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
2168 2583 : return z;
2169 :
2170 70 : case t_QUAD: z = cgetg(4,t_QUAD);
2171 70 : p1 = gel(x,1);
2172 70 : gel(z,1) = ZX_copy(p1); av = avma;
2173 70 : p2 = gsqr(gel(x,2));
2174 70 : p3 = gsqr(gel(x,3));
2175 70 : p4 = gmul(gneg_i(gel(p1,2)),p3);
2176 :
2177 70 : if (gequal0(gel(p1,3)))
2178 : {
2179 7 : tetpil = avma;
2180 7 : gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
2181 7 : av = avma;
2182 7 : p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
2183 7 : gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
2184 : }
2185 :
2186 63 : p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
2187 63 : tetpil = avma;
2188 63 : gel(z,2) = gadd(p2,p4);
2189 63 : gel(z,3) = gadd(p1,p3);
2190 63 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
2191 :
2192 46880 : case t_POLMOD:
2193 46880 : return sqr_polmod(gel(x,1), gel(x,2));
2194 :
2195 2933018 : case t_FFELT: return FF_sqr(x);
2196 :
2197 756897 : case t_POL: return RgX_sqr(x);
2198 :
2199 27552 : case t_SER:
2200 27552 : lx = lg(x);
2201 27552 : if (ser_isexactzero(x)) {
2202 21 : GEN z = gcopy(x);
2203 21 : setvalp(z, 2*valp(x));
2204 21 : return z;
2205 : }
2206 27531 : if (lx < 40)
2207 27272 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2208 : else
2209 : {
2210 259 : pari_sp av = avma;
2211 259 : GEN z = init_ser(lx, varn(x), 2*valp(x));
2212 259 : x = ser2pol_i(x, lx);
2213 259 : x = RgXn_sqr(x, lx-2);
2214 259 : return gerepilecopy(av, fill_ser(z,x));
2215 : }
2216 :
2217 14 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2218 14 : gel(z,1) = gsqr(gel(x,1));
2219 14 : gel(z,2) = gsqr(gel(x,2)); return z;
2220 :
2221 980 : case t_MAT: return RgM_sqr(x);
2222 95617 : case t_VEC: /* handle extended t_QFB */
2223 95617 : case t_QFB: return qfbsqr(x);
2224 658 : case t_VECSMALL:
2225 658 : z = cgetg_copy(x, &lx);
2226 16289 : for (i=1; i<lx; i++)
2227 : {
2228 15631 : long xi = x[i];
2229 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2230 15631 : z[i] = x[xi];
2231 : }
2232 658 : return z;
2233 : }
2234 6 : pari_err_TYPE2("*",x,x);
2235 : return NULL; /* LCOV_EXCL_LINE */
2236 : }
2237 :
2238 : /********************************************************************/
2239 : /** **/
2240 : /** DIVISION **/
2241 : /** **/
2242 : /********************************************************************/
2243 : static GEN
2244 30832 : div_rfrac_scal(GEN x, GEN y)
2245 : {
2246 30832 : pari_sp av = avma;
2247 30832 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2248 30832 : return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
2249 : }
2250 : static GEN
2251 37459 : div_scal_rfrac(GEN x, GEN y)
2252 : {
2253 37459 : GEN y1 = gel(y,1), y2 = gel(y,2);
2254 37459 : pari_sp av = avma;
2255 37459 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2256 : {
2257 7 : if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
2258 0 : y1 = gel(y1,2);
2259 : }
2260 37452 : return RgX_Rg_mul(y2, gdiv(x,y1));
2261 : }
2262 : static GEN
2263 1186635 : div_rfrac(GEN x, GEN y)
2264 1186635 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2265 :
2266 : /* x != 0 */
2267 : static GEN
2268 1290013 : div_ser_scal(GEN y, GEN x)
2269 : {
2270 : long i, l;
2271 : GEN z;
2272 1290013 : if (ser_isexactzero(y))
2273 : {
2274 21 : z = scalarser(lg(y) == 2? Rg_get_0(x): gdiv(gel(y,2), x), varn(y), 1);
2275 21 : setvalp(z, valp(y)); return z;
2276 : }
2277 1289992 : z = cgetg_copy(y, &l); z[1] = y[1];
2278 5913412 : for (i = 2; i < l; i++) gel(z,i) = gdiv(gel(y,i), x);
2279 1289992 : return normalizeser(z);
2280 : }
2281 : GEN
2282 7553 : ser_normalize(GEN x)
2283 : {
2284 7553 : long i, lx = lg(x);
2285 : GEN c, z;
2286 7553 : if (lx == 2) return x;
2287 7553 : c = gel(x,2); if (gequal1(c)) return x;
2288 7476 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2289 107835 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2290 7476 : return z;
2291 : }
2292 :
2293 : /* y != 0 */
2294 : static GEN
2295 4819083 : div_T_scal(GEN x, GEN y, long tx) {
2296 4819083 : switch(tx)
2297 : {
2298 3502112 : case t_POL: return RgX_Rg_div(x, y);
2299 1290006 : case t_SER: return div_ser_scal(x, y);
2300 26968 : case t_RFRAC: return div_rfrac_scal(x,y);
2301 : }
2302 0 : pari_err_TYPE2("/",x,y);
2303 : return NULL; /* LCOV_EXCL_LINE */
2304 : }
2305 :
2306 : static GEN
2307 9257631 : div_scal_pol(GEN x, GEN y) {
2308 9257631 : long ly = lg(y);
2309 : pari_sp av;
2310 9257631 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2311 9179356 : if (isrationalzero(x)) return zeropol(varn(y));
2312 7130116 : av = avma;
2313 7130116 : return gerepileupto(av, gred_rfrac_simple(x,y));
2314 : }
2315 : static GEN
2316 18326 : div_scal_ser(GEN x, GEN y)
2317 18326 : { pari_sp av = avma; return gerepileupto(av, gmul(x, ser_inv(y))); }
2318 : static GEN
2319 9266194 : div_scal_T(GEN x, GEN y, long ty) {
2320 9266194 : switch(ty)
2321 : {
2322 9211172 : case t_POL: return div_scal_pol(x, y);
2323 18319 : case t_SER: return div_scal_ser(x, y);
2324 36703 : case t_RFRAC: return div_scal_rfrac(x, y);
2325 : }
2326 0 : pari_err_TYPE2("/",x,y);
2327 : return NULL; /* LCOV_EXCL_LINE */
2328 : }
2329 :
2330 : /* assume tx = ty = t_SER, same variable vx */
2331 : static GEN
2332 753667 : div_ser(GEN x, GEN y, long vx)
2333 : {
2334 753667 : long e, v = valp(x) - valp(y), lx = lg(x), ly = lg(y);
2335 753667 : GEN y0 = y, z;
2336 753667 : pari_sp av = avma;
2337 :
2338 753667 : if (!signe(y)) pari_err_INV("div_ser", y);
2339 753660 : if (ser_isexactzero(x))
2340 : {
2341 59934 : if (lx == 2) return zeroser(vx, v);
2342 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2343 147 : setvalp(z, v); return z;
2344 : }
2345 693726 : if (lx < ly) ly = lx;
2346 693726 : y = ser2pol_approx(y, ly, &e);
2347 693726 : if (e) { v -= e; ly -= e; if (ly <= 2) pari_err_INV("div_ser", y0); }
2348 693726 : z = init_ser(ly, vx, v);
2349 693726 : if (ly == 3)
2350 : {
2351 15658 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2352 15658 : return gerepileupto(av, z);
2353 : }
2354 678068 : x = ser2pol_i(x, ly);
2355 678068 : y = RgXn_div_i(x, y, ly-2);
2356 678068 : return gerepilecopy(av, fill_ser(z,y));
2357 : }
2358 : /* x,y compatible PADIC */
2359 : static GEN
2360 2814678 : divpp(GEN x, GEN y) {
2361 : pari_sp av;
2362 : long a, b;
2363 : GEN z, M;
2364 :
2365 2814678 : if (!signe(gel(y,4))) pari_err_INV("divpp",y);
2366 2814680 : if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
2367 2814365 : a = precp(x);
2368 2814365 : b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
2369 2814365 : z = cgetg(5, t_PADIC);
2370 2814069 : z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
2371 2813971 : gel(z,2) = icopy(gel(x,2));
2372 2813794 : gel(z,3) = icopy(M); av = avma;
2373 2813885 : gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
2374 2814639 : return z;
2375 : }
2376 : static GEN
2377 27405 : div_polmod_same(GEN T, GEN x, GEN y)
2378 : {
2379 27405 : long v = varn(T);
2380 27405 : GEN a, z = cgetg(3, t_POLMOD);
2381 27405 : gel(z,1) = RgX_copy(T);
2382 27405 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2383 9772 : a = gdiv(x, y);
2384 17633 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2385 119 : {
2386 119 : pari_sp av = avma;
2387 119 : a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
2388 : }
2389 17514 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2390 5145 : {
2391 5145 : pari_sp av = avma;
2392 5145 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2393 5145 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2394 5145 : a = gerepileupto(av, a);
2395 : }
2396 : else
2397 : {
2398 12369 : pari_sp av = avma;
2399 12369 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2400 12369 : a = gerepileupto(av, a);
2401 : }
2402 27405 : gel(z,2) = a; return z;
2403 : }
2404 : GEN
2405 374042008 : gdiv(GEN x, GEN y)
2406 : {
2407 374042008 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2408 : pari_sp av, tetpil;
2409 : GEN z, p1, p2;
2410 :
2411 374042008 : if (tx == ty) switch(tx)
2412 : {
2413 78044212 : case t_INT:
2414 78044212 : return Qdivii(x,y);
2415 :
2416 81952067 : case t_REAL: return divrr(x,y);
2417 21083 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2418 21083 : z = cgetg(3,t_INTMOD);
2419 21083 : if (X==Y || equalii(X,Y))
2420 21076 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2421 7 : gel(z,1) = gcdii(X,Y);
2422 7 : warn_coercion(X,Y,gel(z,1));
2423 7 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2424 7 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
2425 : }
2426 3683729 : case t_FRAC: {
2427 3683729 : GEN x1 = gel(x,1), x2 = gel(x,2);
2428 3683729 : GEN y1 = gel(y,1), y2 = gel(y,2);
2429 3683729 : z = cgetg(3, t_FRAC);
2430 3683729 : p1 = gcdii(x1, y1);
2431 3683728 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2432 3683728 : p1 = gcdii(x2, y2);
2433 3683728 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2434 3683724 : tetpil = avma;
2435 3683724 : gel(z,2) = mulii(x2,y1);
2436 3683723 : gel(z,1) = mulii(x1,y2);
2437 3683726 : normalize_frac(z);
2438 3683727 : fix_frac_if_int_GC(z,tetpil);
2439 3683729 : return z;
2440 : }
2441 7697710 : case t_COMPLEX:
2442 7697710 : if (isintzero(gel(y,1)))
2443 : {
2444 110362 : y = gel(y,2);
2445 110362 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2446 93359 : z = cgetg(3,t_COMPLEX);
2447 93359 : gel(z,1) = gdiv(gel(x,2), y);
2448 93359 : av = avma;
2449 93359 : gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
2450 93359 : return z;
2451 : }
2452 7587346 : av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
2453 7587414 : return gerepile(av, tetpil, gdiv(p2,p1));
2454 :
2455 767638 : case t_PADIC:
2456 767638 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
2457 767638 : return divpp(x, y);
2458 :
2459 238 : case t_QUAD:
2460 238 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2461 224 : av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
2462 224 : return gerepile(av, tetpil, gdiv(p2,p1));
2463 :
2464 236484 : case t_FFELT: return FF_div(x,y);
2465 :
2466 31780 : case t_POLMOD:
2467 31780 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2468 27405 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2469 : else {
2470 4375 : av = avma;
2471 4375 : z = gerepileupto(av, gmul(x, ginv(y)));
2472 : }
2473 31780 : return z;
2474 :
2475 18459397 : case t_POL:
2476 18459397 : vx = varn(x);
2477 18459397 : vy = varn(y);
2478 18459397 : if (vx != vy) {
2479 102557 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2480 46459 : else return div_scal_pol(x, y);
2481 : }
2482 18356840 : if (!signe(y)) pari_err_INV("gdiv",y);
2483 18356840 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2484 18233428 : av = avma;
2485 18233428 : return gerepileupto(av, gred_rfrac2(x,y));
2486 :
2487 753688 : case t_SER:
2488 753688 : vx = varn(x);
2489 753688 : vy = varn(y);
2490 753688 : if (vx != vy) {
2491 21 : if (varncmp(vx, vy) < 0)
2492 : {
2493 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2494 7 : return div_ser_scal(x, y);
2495 : }
2496 7 : return div_scal_ser(x, y);
2497 : }
2498 753667 : return div_ser(x, y, vx);
2499 1191115 : case t_RFRAC:
2500 1191115 : vx = varn(gel(x,2));
2501 1191115 : vy = varn(gel(y,2));
2502 1191115 : if (vx != vy) {
2503 4480 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2504 756 : else return div_scal_rfrac(x, y);
2505 : }
2506 1186635 : return div_rfrac(x,y);
2507 :
2508 21 : case t_VEC: /* handle extended t_QFB */
2509 21 : case t_QFB: av = avma; return gerepileupto(av, qfbcomp(x, ginv(y)));
2510 :
2511 14 : case t_MAT:
2512 14 : av = avma; p1 = RgM_inv(y);
2513 14 : if (!p1) pari_err_INV("gdiv",y);
2514 14 : return gerepileupto(av, RgM_mul(x, p1));
2515 :
2516 0 : default: pari_err_TYPE2("/",x,y);
2517 : }
2518 :
2519 181208208 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2520 : {
2521 3845062 : long s = signe(x);
2522 3845062 : if (!s) {
2523 476184 : if (gequal0(y)) pari_err_INV("gdiv",y);
2524 476184 : switch (ty)
2525 : {
2526 473027 : default: return gen_0;
2527 21 : case t_INTMOD:
2528 21 : z = cgetg(3,t_INTMOD);
2529 21 : gel(z,1) = icopy(gel(y,1));
2530 21 : gel(z,2) = gen_0; return z;
2531 3136 : case t_FFELT: return FF_zero(y);
2532 : }
2533 : }
2534 3368878 : if (is_pm1(x)) {
2535 1095008 : if (s > 0) return ginv(y);
2536 189782 : av = avma; return gerepileupto(av, ginv(gneg(y)));
2537 : }
2538 2273870 : switch(ty)
2539 : {
2540 798458 : case t_REAL: return divir(x,y);
2541 21 : case t_INTMOD:
2542 21 : z = cgetg(3, t_INTMOD);
2543 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2544 780681 : case t_FRAC:
2545 780681 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2546 780681 : if (equali1(p1))
2547 : {
2548 249857 : set_avma((pari_sp)z);
2549 249857 : gel(z,2) = icopy(gel(y,1));
2550 249857 : gel(z,1) = mulii(gel(y,2), x);
2551 249857 : normalize_frac(z);
2552 249857 : fix_frac_if_int(z);
2553 : }
2554 : else
2555 : {
2556 530824 : x = diviiexact(x,p1); tetpil = avma;
2557 530824 : gel(z,2) = diviiexact(gel(y,1), p1);
2558 530824 : gel(z,1) = mulii(gel(y,2), x);
2559 530824 : normalize_frac(z);
2560 530824 : fix_frac_if_int_GC(z,tetpil);
2561 : }
2562 780681 : return z;
2563 :
2564 273 : case t_FFELT: return Z_FF_div(x,y);
2565 692713 : case t_COMPLEX: return divRc(x,y);
2566 1505 : case t_PADIC: return divTp(x, y);
2567 231 : case t_QUAD:
2568 231 : av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
2569 231 : return gerepile(av, tetpil, gdiv(p2,p1));
2570 : }
2571 : }
2572 177363133 : if (gequal0(y))
2573 : {
2574 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2575 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2576 : }
2577 :
2578 177373852 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2579 : {
2580 142413540 : case t_REAL:
2581 142413540 : switch(ty)
2582 : {
2583 140087005 : case t_INT: return divri(x,y);
2584 507168 : case t_FRAC:
2585 507168 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2586 507168 : return gerepileuptoleaf(av, z);
2587 1819430 : case t_COMPLEX: return divRc(x, y);
2588 42 : case t_QUAD: return divfq(x, y, realprec(x));
2589 18 : default: pari_err_TYPE2("/",x,y);
2590 : }
2591 :
2592 287 : case t_INTMOD:
2593 : switch(ty)
2594 : {
2595 203 : case t_INT:
2596 203 : z = cgetg(3, t_INTMOD);
2597 203 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2598 28 : case t_FRAC: { GEN X = gel(x,1);
2599 28 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2600 28 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2601 : }
2602 7 : case t_FFELT:
2603 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2604 0 : pari_err_OP("/",x,y);
2605 7 : return Z_FF_div(gel(x,2),y);
2606 :
2607 0 : case t_COMPLEX:
2608 0 : av = avma;
2609 0 : return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
2610 :
2611 14 : case t_QUAD:
2612 14 : av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
2613 14 : return gerepile(av,tetpil, gdiv(p2,p1));
2614 :
2615 7 : case t_PADIC: { GEN X = gel(x,1);
2616 7 : z = cgetg(3, t_INTMOD);
2617 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2618 : }
2619 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2620 : }
2621 :
2622 : case t_FRAC:
2623 : switch(ty)
2624 : {
2625 1998009 : case t_INT: z = cgetg(3, t_FRAC);
2626 1998009 : p1 = gcdii(y,gel(x,1));
2627 1998009 : if (equali1(p1))
2628 : {
2629 755923 : set_avma((pari_sp)z); tetpil = 0;
2630 755923 : gel(z,1) = icopy(gel(x,1));
2631 : }
2632 : else
2633 : {
2634 1242086 : y = diviiexact(y,p1); tetpil = avma;
2635 1242086 : gel(z,1) = diviiexact(gel(x,1), p1);
2636 : }
2637 1998009 : gel(z,2) = mulii(gel(x,2),y);
2638 1998009 : normalize_frac(z);
2639 1998009 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2640 1998009 : return z;
2641 :
2642 57745 : case t_REAL:
2643 57745 : av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
2644 57745 : return gerepile(av, tetpil, divir(gel(x,1), p1));
2645 :
2646 7 : case t_INTMOD: { GEN Y = gel(y,1);
2647 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2648 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2649 : }
2650 :
2651 28 : case t_FFELT: av=avma;
2652 28 : return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2653 :
2654 483 : case t_COMPLEX: return divRc(x, y);
2655 :
2656 2128 : case t_PADIC:
2657 2128 : if (!signe(gel(x,1))) return gen_0;
2658 2128 : return divTp(x, y);
2659 :
2660 42 : case t_QUAD:
2661 42 : av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
2662 42 : return gerepile(av,tetpil,gdiv(p2,p1));
2663 : }
2664 :
2665 : case t_FFELT:
2666 133 : switch (ty)
2667 : {
2668 49 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2669 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2670 7 : case t_INTMOD:
2671 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2672 0 : pari_err_OP("/",x,y);
2673 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2674 49 : default:
2675 49 : pari_err_TYPE2("/",x,y);
2676 : }
2677 0 : break;
2678 :
2679 11942807 : case t_COMPLEX:
2680 : switch(ty)
2681 : {
2682 11942788 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2683 0 : case t_INTMOD: return mulRc_direct(ginv(y), x);
2684 0 : case t_PADIC:
2685 0 : return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
2686 0 : case t_QUAD:
2687 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2688 0 : return divfq(x, y, lx);
2689 : }
2690 :
2691 : case t_PADIC:
2692 : switch(ty)
2693 : {
2694 1327291 : case t_INT: case t_FRAC: { GEN p = gel(x,2);
2695 1327200 : return signe(gel(x,4))? divpT(x, y)
2696 2654491 : : zeropadic(p, valp(x) - Q_pval(y,p));
2697 : }
2698 7 : case t_INTMOD: { GEN Y = gel(y,1);
2699 7 : z = cgetg(3, t_INTMOD);
2700 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2701 : }
2702 14 : case t_COMPLEX: case t_QUAD:
2703 14 : av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
2704 14 : return gerepile(av,tetpil,gdiv(p1,p2));
2705 :
2706 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2707 : }
2708 :
2709 : case t_QUAD:
2710 : switch (ty)
2711 : {
2712 1190 : case t_INT: case t_INTMOD: case t_FRAC:
2713 1190 : z = cgetg(4,t_QUAD);
2714 1190 : gel(z,1) = ZX_copy(gel(x,1));
2715 1190 : gel(z,2) = gdiv(gel(x,2), y);
2716 1190 : gel(z,3) = gdiv(gel(x,3), y); return z;
2717 63 : case t_REAL: return divqf(x, y, realprec(y));
2718 14 : case t_PADIC: return divTp(x, y);
2719 0 : case t_COMPLEX:
2720 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2721 0 : return divqf(x, y, ly);
2722 : }
2723 : }
2724 19629865 : switch(ty) {
2725 374309 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2726 374309 : return gmul(x, ginv(y)); /* missing gerepile, for speed */
2727 42 : case t_MAT:
2728 42 : av = avma; p1 = RgM_inv(y);
2729 35 : if (!p1) pari_err_INV("gdiv",y);
2730 28 : return gerepileupto(av, gmul(x, p1));
2731 0 : case t_VEC: case t_COL:
2732 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2733 0 : pari_err_TYPE2("/",x,y);
2734 : }
2735 19255514 : switch(tx) {
2736 3169904 : case t_VEC: case t_COL: case t_MAT:
2737 3169904 : z = cgetg_copy(x, &lx);
2738 11081150 : for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
2739 3169798 : return z;
2740 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2741 0 : pari_err_TYPE2("/",x,y);
2742 : }
2743 :
2744 16085610 : vy = gvar(y);
2745 16086422 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2746 207456 : vx = varn(X);
2747 207456 : if (vx != vy) {
2748 206882 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2749 206553 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2750 : }
2751 : /* y is POL, SER or RFRAC */
2752 574 : av = avma;
2753 574 : switch(ty)
2754 : {
2755 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2756 574 : default: y = ginvmod(gmod(y,X), X);
2757 : }
2758 567 : return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
2759 : }
2760 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2761 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2762 : * variable NO_VARIABLE, then the operation is incorrect */
2763 15878966 : vx = gvar(x);
2764 15878961 : if (vx != vy) { /* includes cases where one is scalar */
2765 14084948 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2766 9265863 : else return div_scal_T(x, y, ty);
2767 : }
2768 1794013 : switch(tx)
2769 : {
2770 1046143 : case t_POL:
2771 : switch(ty)
2772 : {
2773 28 : case t_SER:
2774 : {
2775 28 : GEN y0 = y;
2776 : long v;
2777 28 : ly = lg(y); /* > 2 */
2778 28 : av = avma; v = RgX_valrem(x, &x);
2779 28 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2780 14 : v -= valp(y);
2781 14 : y = ser2pol_approx(y, ly, &i); if (i) { ly -= i; v -= i; }
2782 14 : if (ly == 2) pari_err_INV("gdiv", y0);
2783 7 : z = init_ser(ly, vx, v);
2784 7 : return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, ly-2)));
2785 : }
2786 :
2787 1046115 : case t_RFRAC:
2788 : {
2789 1046115 : GEN y1 = gel(y,1), y2 = gel(y,2);
2790 1046115 : if (typ(y1) == t_POL && varn(y1) == vx)
2791 35 : return mul_rfrac_scal(y2, y1, x);
2792 1046080 : av = avma;
2793 1046080 : return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2794 : }
2795 : }
2796 0 : break;
2797 :
2798 478161 : case t_SER:
2799 : switch(ty)
2800 : {
2801 478147 : case t_POL:
2802 : {
2803 478147 : long v = valp(x);
2804 478147 : lx = lg(x);
2805 478147 : if (lx == 2) return zeroser(vx, v - RgX_val(y));
2806 478042 : av = avma;
2807 478042 : x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
2808 478042 : z = init_ser(lx, vx, v);
2809 478042 : if (!signe(x)) setsigne(z,0);
2810 478042 : return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
2811 : }
2812 14 : case t_RFRAC:
2813 14 : av = avma;
2814 14 : return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2815 : }
2816 0 : break;
2817 :
2818 269687 : case t_RFRAC:
2819 : switch(ty)
2820 : {
2821 269687 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2822 0 : case t_SER:
2823 0 : av = avma;
2824 0 : return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2825 : }
2826 0 : break;
2827 : }
2828 22 : pari_err_TYPE2("/",x,y);
2829 : return NULL; /* LCOV_EXCL_LINE */
2830 : }
2831 :
2832 : /********************************************************************/
2833 : /** **/
2834 : /** SIMPLE MULTIPLICATION **/
2835 : /** **/
2836 : /********************************************************************/
2837 : GEN
2838 246693993 : gmulsg(long s, GEN y)
2839 : {
2840 : long ly, i;
2841 : pari_sp av;
2842 : GEN z;
2843 :
2844 246693993 : switch(typ(y))
2845 : {
2846 156937086 : case t_INT: return mulsi(s,y);
2847 69151131 : case t_REAL: return s? mulsr(s,y): gen_0; /* gmul semantic */
2848 174905 : case t_INTMOD: { GEN p = gel(y,1);
2849 174905 : z = cgetg(3,t_INTMOD);
2850 174907 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
2851 174907 : gel(z,1) = icopy(p); return z;
2852 : }
2853 547745 : case t_FFELT: return FF_Z_mul(y,stoi(s));
2854 5292967 : case t_FRAC:
2855 5292967 : if (!s) return gen_0;
2856 5288053 : z = cgetg(3,t_FRAC);
2857 5288053 : i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
2858 5288053 : if (i == 1)
2859 : {
2860 2021662 : gel(z,2) = icopy(gel(y,2));
2861 2021662 : gel(z,1) = mulis(gel(y,1), s);
2862 : }
2863 : else
2864 : {
2865 3266391 : gel(z,2) = diviuexact(gel(y,2), (ulong)i);
2866 3266391 : gel(z,1) = mulis(gel(y,1), s/i);
2867 3266391 : fix_frac_if_int(z);
2868 : }
2869 5288053 : return z;
2870 :
2871 13102396 : case t_COMPLEX:
2872 13102396 : if (!s) return gen_0;
2873 13028048 : z = cgetg(3, t_COMPLEX);
2874 13027861 : gel(z,1) = gmulsg(s,gel(y,1));
2875 13025930 : gel(z,2) = gmulsg(s,gel(y,2)); return z;
2876 :
2877 1435 : case t_PADIC:
2878 1435 : if (!s) return gen_0;
2879 1435 : av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
2880 :
2881 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2882 7 : gel(z,1) = ZX_copy(gel(y,1));
2883 7 : gel(z,2) = gmulsg(s,gel(y,2));
2884 7 : gel(z,3) = gmulsg(s,gel(y,3)); return z;
2885 :
2886 74802 : case t_POLMOD:
2887 74802 : retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
2888 :
2889 585875 : case t_POL:
2890 585875 : if (!signe(y)) return RgX_copy(y);
2891 572295 : if (!s) return scalarpol(Rg_get_0(y), varn(y));
2892 570671 : z = cgetg_copy(y, &ly); z[1]=y[1];
2893 2342076 : for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2894 570671 : return normalizepol_lg(z, ly);
2895 :
2896 182 : case t_SER:
2897 182 : if (ser_isexactzero(y)) return gcopy(y);
2898 182 : if (!s) return Rg_get_0(y);
2899 182 : z = cgetg_copy(y, &ly); z[1]=y[1];
2900 3864 : for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2901 182 : return normalizeser(z);
2902 :
2903 14 : case t_RFRAC:
2904 14 : if (!s) return zeropol(varn(gel(y,2)));
2905 14 : if (s == 1) return gcopy(y);
2906 14 : if (s == -1) return gneg(y);
2907 14 : return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
2908 :
2909 837166 : case t_VEC: case t_COL: case t_MAT:
2910 837166 : z = cgetg_copy(y, &ly);
2911 2480754 : for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2912 837154 : return z;
2913 : }
2914 0 : pari_err_TYPE("gmulsg",y);
2915 : return NULL; /* LCOV_EXCL_LINE */
2916 : }
2917 :
2918 : GEN
2919 132364991 : gmulug(ulong s, GEN y)
2920 : {
2921 : long ly, i;
2922 : pari_sp av;
2923 : GEN z;
2924 :
2925 132364991 : switch(typ(y))
2926 : {
2927 130456383 : case t_INT: return mului(s,y);
2928 1046319 : case t_REAL: return s? mulur(s,y): gen_0; /* gmul semantic */
2929 364 : case t_INTMOD: { GEN p = gel(y,1);
2930 364 : z = cgetg(3,t_INTMOD);
2931 364 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(mului(s,gel(y,2)), p));
2932 364 : gel(z,1) = icopy(p); return z;
2933 : }
2934 413 : case t_FFELT: return FF_Z_mul(y,utoi(s));
2935 802877 : case t_FRAC:
2936 802877 : if (!s) return gen_0;
2937 802863 : z = cgetg(3,t_FRAC);
2938 802863 : i = ugcd(s, umodiu(gel(y,2), s));
2939 802863 : if (i == 1)
2940 : {
2941 665972 : gel(z,2) = icopy(gel(y,2));
2942 665972 : gel(z,1) = muliu(gel(y,1), s);
2943 : }
2944 : else
2945 : {
2946 136891 : gel(z,2) = diviuexact(gel(y,2), i);
2947 136891 : gel(z,1) = muliu(gel(y,1), s/i);
2948 136891 : fix_frac_if_int(z);
2949 : }
2950 802863 : return z;
2951 :
2952 24750 : case t_COMPLEX:
2953 24750 : if (!s) return gen_0;
2954 24750 : z = cgetg(3, t_COMPLEX);
2955 24750 : gel(z,1) = gmulug(s,gel(y,1));
2956 24750 : gel(z,2) = gmulug(s,gel(y,2)); return z;
2957 :
2958 3898 : case t_PADIC:
2959 3898 : if (!s) return gen_0;
2960 3898 : av = avma; return gerepileupto(av, mulpp(cvtop2(utoi(s),y), y));
2961 :
2962 0 : case t_QUAD: z = cgetg(4, t_QUAD);
2963 0 : gel(z,1) = ZX_copy(gel(y,1));
2964 0 : gel(z,2) = gmulug(s,gel(y,2));
2965 0 : gel(z,3) = gmulug(s,gel(y,3)); return z;
2966 :
2967 3199 : case t_POLMOD:
2968 3199 : retmkpolmod(gmulug(s,gel(y,2)), RgX_copy(gel(y,1)));
2969 :
2970 12642 : case t_POL:
2971 12642 : if (!signe(y)) return RgX_copy(y);
2972 11865 : if (!s) return scalarpol(Rg_get_0(y), varn(y));
2973 11865 : z = cgetg_copy(y, &ly); z[1]=y[1];
2974 34923 : for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
2975 11865 : return normalizepol_lg(z, ly);
2976 :
2977 0 : case t_SER:
2978 0 : if (ser_isexactzero(y)) return gcopy(y);
2979 0 : if (!s) return Rg_get_0(y);
2980 0 : z = cgetg_copy(y, &ly); z[1]=y[1];
2981 0 : for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
2982 0 : return normalizeser(z);
2983 :
2984 0 : case t_RFRAC:
2985 0 : if (!s) return zeropol(varn(gel(y,2)));
2986 0 : if (s == 1) return gcopy(y);
2987 0 : return mul_rfrac_scal(gel(y,1), gel(y,2), utoi(s));
2988 :
2989 14147 : case t_VEC: case t_COL: case t_MAT:
2990 14147 : z = cgetg_copy(y, &ly);
2991 73922 : for (i=1; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
2992 14147 : return z;
2993 : }
2994 0 : pari_err_TYPE("gmulsg",y);
2995 : return NULL; /* LCOV_EXCL_LINE */
2996 : }
2997 :
2998 : /********************************************************************/
2999 : /** **/
3000 : /** SIMPLE DIVISION **/
3001 : /** **/
3002 : /********************************************************************/
3003 :
3004 : GEN
3005 7219500 : gdivgs(GEN x, long s)
3006 : {
3007 7219500 : long tx = typ(x), lx, i;
3008 : pari_sp av;
3009 : GEN z;
3010 :
3011 7219500 : if (!s)
3012 : {
3013 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3014 0 : pari_err_INV("gdivgs",gen_0);
3015 : }
3016 7219517 : switch(tx)
3017 : {
3018 1326121 : case t_INT: return Qdivis(x, s);
3019 4690896 : case t_REAL: return divrs(x,s);
3020 :
3021 357 : case t_INTMOD:
3022 357 : z = cgetg(3, t_INTMOD);
3023 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3024 :
3025 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3026 :
3027 383917 : case t_FRAC: z = cgetg(3, t_FRAC);
3028 383917 : i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
3029 383917 : if (i == 1)
3030 : {
3031 294265 : gel(z,2) = mulsi(s, gel(x,2));
3032 294265 : gel(z,1) = icopy(gel(x,1));
3033 : }
3034 : else
3035 : {
3036 89652 : gel(z,2) = mulsi(s/i, gel(x,2));
3037 89652 : gel(z,1) = divis(gel(x,1), i);
3038 : }
3039 383917 : normalize_frac(z);
3040 383917 : fix_frac_if_int(z); return z;
3041 :
3042 756423 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3043 756423 : gel(z,1) = gdivgs(gel(x,1),s);
3044 756421 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3045 :
3046 133 : case t_PADIC: /* divpT */
3047 : {
3048 133 : GEN p = gel(x,2);
3049 133 : if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
3050 133 : av = avma;
3051 133 : return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
3052 : }
3053 :
3054 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3055 28 : gel(z,1) = ZX_copy(gel(x,1));
3056 28 : gel(z,2) = gdivgs(gel(x,2),s);
3057 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3058 :
3059 29659 : case t_POLMOD:
3060 29659 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3061 :
3062 91 : case t_RFRAC:
3063 91 : if (s == 1) return gcopy(x);
3064 84 : else if (s == -1) return gneg(x);
3065 84 : return div_rfrac_scal(x, stoi(s));
3066 :
3067 29442 : case t_POL: case t_SER:
3068 29442 : z = cgetg_copy(x, &lx); z[1] = x[1];
3069 117191 : for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
3070 29442 : return z;
3071 1715 : case t_VEC: case t_COL: case t_MAT:
3072 1715 : z = cgetg_copy(x, &lx);
3073 40971 : for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
3074 1715 : return z;
3075 :
3076 : }
3077 0 : pari_err_TYPE2("/",x, stoi(s));
3078 : return NULL; /* LCOV_EXCL_LINE */
3079 : }
3080 :
3081 : GEN
3082 42784391 : gdivgu(GEN x, ulong s)
3083 : {
3084 42784391 : long tx = typ(x), lx, i;
3085 : pari_sp av;
3086 : GEN z;
3087 :
3088 42784391 : if (!s)
3089 : {
3090 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3091 0 : pari_err_INV("gdivgu",gen_0);
3092 : }
3093 42784327 : switch(tx)
3094 : {
3095 16943589 : case t_INT: return Qdiviu(x, s);
3096 12889670 : case t_REAL: return divru(x,s);
3097 :
3098 210315 : case t_INTMOD:
3099 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3100 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3101 :
3102 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3103 :
3104 280325 : case t_FRAC: z = cgetg(3, t_FRAC);
3105 280325 : i = ugcd(s, umodiu(gel(x,1), s));
3106 280325 : if (i == 1)
3107 : {
3108 148028 : gel(z,2) = mului(s, gel(x,2));
3109 148028 : gel(z,1) = icopy(gel(x,1));
3110 : }
3111 : else
3112 : {
3113 132297 : gel(z,2) = mului(s/i, gel(x,2));
3114 132297 : gel(z,1) = divis(gel(x,1), i);
3115 : }
3116 280325 : normalize_frac(z);
3117 280325 : fix_frac_if_int(z); return z;
3118 :
3119 11556147 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3120 11556148 : gel(z,1) = gdivgu(gel(x,1),s);
3121 11556146 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3122 :
3123 850395 : case t_PADIC: /* divpT */
3124 : {
3125 850395 : GEN p = gel(x,2);
3126 850395 : if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
3127 716264 : av = avma;
3128 716264 : return gerepileupto(av, divpp(x, cvtop2(utoi(s),x)));
3129 : }
3130 :
3131 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3132 0 : gel(z,1) = ZX_copy(gel(x,1));
3133 0 : gel(z,2) = gdivgu(gel(x,2),s);
3134 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3135 :
3136 1456 : case t_POLMOD:
3137 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3138 :
3139 56 : case t_RFRAC:
3140 56 : if (s == 1) return gcopy(x);
3141 56 : return div_rfrac_scal(x, utoi(s));
3142 :
3143 51765 : case t_POL: case t_SER:
3144 51765 : z = cgetg_copy(x, &lx); z[1] = x[1];
3145 228326 : for (i=2; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
3146 51765 : return z;
3147 301 : case t_VEC: case t_COL: case t_MAT:
3148 301 : z = cgetg_copy(x, &lx);
3149 1092 : for (i=1; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
3150 301 : return z;
3151 : }
3152 0 : pari_err_TYPE2("/",x, utoi(s));
3153 : return NULL; /* LCOV_EXCL_LINE */
3154 : }
3155 :
3156 : /* x / (i*(i+1)) */
3157 : GEN
3158 105279624 : divrunextu(GEN x, ulong i)
3159 : {
3160 105279624 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3161 0 : return divri(x, muluu(i , i+1));
3162 : else
3163 105279624 : return divru(x, i*(i+1));
3164 : }
3165 : /* x / (i*(i+1)) */
3166 : GEN
3167 804801 : gdivgunextu(GEN x, ulong i)
3168 : {
3169 804801 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3170 0 : return gdivgu(x, i*(i+1));
3171 : else
3172 804801 : return gdiv(x, muluu(i, i+1));
3173 : }
3174 :
3175 : /* True shift (exact multiplication by 2^n) */
3176 : GEN
3177 85971815 : gmul2n(GEN x, long n)
3178 : {
3179 : long lx, i, k, l;
3180 : GEN z, a, b;
3181 :
3182 85971815 : switch(typ(x))
3183 : {
3184 19111446 : case t_INT:
3185 19111446 : if (n>=0) return shifti(x,n);
3186 3383056 : if (!signe(x)) return gen_0;
3187 3291636 : l = vali(x); n = -n;
3188 3291721 : if (n<=l) return shifti(x,-n);
3189 375044 : z = cgetg(3,t_FRAC);
3190 375044 : gel(z,1) = shifti(x,-l);
3191 375044 : gel(z,2) = int2n(n-l); return z;
3192 :
3193 49073329 : case t_REAL:
3194 49073329 : return shiftr(x,n);
3195 :
3196 180800 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3197 180800 : z = cgetg(3,t_INTMOD);
3198 180800 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3199 82639 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
3200 82642 : gel(z,1) = icopy(b); return z;
3201 :
3202 217234 : case t_FFELT: return FF_mul2n(x,n);
3203 :
3204 839637 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3205 839637 : l = vali(a);
3206 839637 : k = vali(b);
3207 839637 : if (n+l >= k)
3208 : {
3209 299999 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3210 236641 : l = n-k; k = -k;
3211 : }
3212 : else
3213 : {
3214 539638 : k = -(l+n); l = -l;
3215 : }
3216 776279 : z = cgetg(3,t_FRAC);
3217 776279 : gel(z,1) = shifti(a,l);
3218 776279 : gel(z,2) = shifti(b,k); return z;
3219 :
3220 12132834 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3221 12132833 : gel(z,1) = gmul2n(gel(x,1),n);
3222 12132836 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3223 :
3224 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3225 105 : gel(z,1) = ZX_copy(gel(x,1));
3226 105 : gel(z,2) = gmul2n(gel(x,2),n);
3227 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3228 :
3229 44030 : case t_POLMOD:
3230 44030 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3231 :
3232 1578729 : case t_POL:
3233 1578729 : z = cgetg_copy(x, &lx); z[1] = x[1];
3234 8586714 : for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3235 1578600 : return normalizepol_lg(z, lx); /* needed if char = 2 */
3236 106158 : case t_SER:
3237 106158 : if (ser_isexactzero(x)) return gcopy(x);
3238 106130 : z = cgetg_copy(x, &lx); z[1] = x[1];
3239 642327 : for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3240 106130 : return normalizeser(z); /* needed if char = 2 */
3241 2684422 : case t_VEC: case t_COL: case t_MAT:
3242 2684422 : z = cgetg_copy(x, &lx);
3243 12155509 : for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3244 2684316 : return z;
3245 :
3246 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3247 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3248 :
3249 4536 : case t_PADIC: /* int2n wrong if n < 0 */
3250 4536 : return gmul(gmul2n(gen_1,n),x);
3251 : }
3252 0 : pari_err_TYPE("gmul2n",x);
3253 : return NULL; /* LCOV_EXCL_LINE */
3254 : }
3255 :
3256 : /*******************************************************************/
3257 : /* */
3258 : /* INVERSE */
3259 : /* */
3260 : /*******************************************************************/
3261 : static GEN
3262 79175 : inv_polmod(GEN T, GEN x)
3263 : {
3264 79175 : GEN z = cgetg(3,t_POLMOD), a;
3265 79173 : gel(z,1) = RgX_copy(T);
3266 79167 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3267 62251 : a = ginv(x);
3268 : else
3269 : {
3270 16916 : if (lg(T) == 5) /* quadratic fields */
3271 12705 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3272 : else
3273 4211 : a = RgXQ_inv(x, T);
3274 : }
3275 79167 : gel(z,2) = a; return z;
3276 : }
3277 :
3278 : GEN
3279 21599985 : ginv(GEN x)
3280 : {
3281 : long s;
3282 : pari_sp av, tetpil;
3283 : GEN z, y, p1, p2;
3284 :
3285 21599985 : switch(typ(x))
3286 : {
3287 3360888 : case t_INT:
3288 3360888 : if (is_pm1(x)) return icopy(x);
3289 2528513 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3290 2528499 : z = cgetg(3,t_FRAC);
3291 2528534 : gel(z,1) = s<0? gen_m1: gen_1;
3292 2528534 : gel(z,2) = absi(x); return z;
3293 :
3294 5265353 : case t_REAL: return invr(x);
3295 :
3296 4641 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3297 4641 : gel(z,1) = icopy(gel(x,1));
3298 4641 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3299 :
3300 386251 : case t_FRAC: {
3301 386251 : GEN a = gel(x,1), b = gel(x,2);
3302 386251 : s = signe(a);
3303 386251 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3304 205765 : z = cgetg(3,t_FRAC);
3305 205765 : gel(z,1) = icopy(b);
3306 205765 : gel(z,2) = icopy(a);
3307 205765 : normalize_frac(z); return z;
3308 : }
3309 6666599 : case t_COMPLEX:
3310 6666599 : av=avma;
3311 6666599 : p1=cxnorm(x);
3312 6665764 : p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
3313 6666475 : tetpil=avma;
3314 6666475 : return gerepile(av,tetpil,divcR(p2,p1));
3315 :
3316 273 : case t_QUAD:
3317 273 : av=avma; p1=quadnorm(x); p2=conj_i(x); tetpil=avma;
3318 273 : return gerepile(av,tetpil,gdiv(p2,p1));
3319 :
3320 346213 : case t_PADIC: z = cgetg(5,t_PADIC);
3321 346213 : if (!signe(gel(x,4))) pari_err_INV("ginv",x);
3322 346206 : z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
3323 346206 : gel(z,2) = icopy(gel(x,2));
3324 346206 : gel(z,3) = icopy(gel(x,3));
3325 346206 : gel(z,4) = Zp_inv(gel(x,4),gel(z,2),precp(x)); return z;
3326 :
3327 79176 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3328 14269 : case t_FFELT: return FF_inv(x);
3329 5321021 : case t_POL: return gred_rfrac_simple(gen_1,x);
3330 18858 : case t_SER: return ser_inv(x);
3331 2926 : case t_RFRAC:
3332 : {
3333 2926 : GEN n = gel(x,1), d = gel(x,2);
3334 2926 : pari_sp av = avma, ltop;
3335 2926 : if (gequal0(n)) pari_err_INV("ginv",x);
3336 :
3337 2926 : n = simplify_shallow(n);
3338 2926 : if (typ(n) != t_POL || varn(n) != varn(d))
3339 : {
3340 2926 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3341 679 : ltop = avma;
3342 679 : z = RgX_Rg_div(d,n);
3343 : } else {
3344 0 : ltop = avma;
3345 0 : z = cgetg(3,t_RFRAC);
3346 0 : gel(z,1) = RgX_copy(d);
3347 0 : gel(z,2) = RgX_copy(n);
3348 : }
3349 679 : stackdummy(av, ltop);
3350 679 : return z;
3351 : }
3352 :
3353 99119 : case t_VEC: /* handle extended t_QFB */
3354 : case t_QFB:
3355 99119 : return qfbpow(x, gen_m1);
3356 34496 : case t_MAT:
3357 34496 : y = RgM_inv(x);
3358 34489 : if (!y) pari_err_INV("ginv",x);
3359 34419 : return y;
3360 28 : case t_VECSMALL:
3361 : {
3362 28 : long i, lx = lg(x)-1;
3363 28 : y = zero_zv(lx);
3364 112 : for (i=1; i<=lx; i++)
3365 : {
3366 84 : long xi = x[i];
3367 84 : if (xi<1 || xi>lx || y[xi])
3368 0 : pari_err_TYPE("ginv [not a permutation]", x);
3369 84 : y[xi] = i;
3370 : }
3371 28 : return y;
3372 : }
3373 : }
3374 0 : pari_err_TYPE("inverse",x);
3375 : return NULL; /* LCOV_EXCL_LINE */
3376 : }
|
