Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** GENERIC OPERATIONS **/
18 : /** (first part) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_mod
25 :
26 : /* assume z[1] was created last */
27 : #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
28 : z = gerepileupto((pari_sp)(z+3), gel(z,1));
29 :
30 : /* assume z[1] was created last */
31 : #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
32 : z = gerepileupto((pari_sp)(z+3), gel(z,1));\
33 : else\
34 : gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
35 :
36 : static void
37 105 : warn_coercion(GEN x, GEN y, GEN z)
38 : {
39 105 : if (DEBUGLEVEL)
40 56 : pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
41 105 : }
42 :
43 : static long
44 28 : kro_quad(GEN x, GEN y)
45 28 : { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
46 :
47 : /* is -1 not a square in Zp, assume p prime */
48 : INLINE int
49 42 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
50 :
51 : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
52 : static GEN mulpp(GEN x, GEN y);
53 : static GEN divpp(GEN x, GEN y);
54 : /* Argument codes for inline routines
55 : * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
56 : * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
57 : * T: some type (to be converted to PADIC)
58 : */
59 : static GEN
60 303695316 : addRc(GEN x, GEN y) {
61 303695316 : GEN z = cgetg(3,t_COMPLEX);
62 303679148 : gel(z,1) = gadd(x,gel(y,1));
63 303687926 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 350776092 : mulRc(GEN x, GEN y) {
67 350776092 : GEN z = cgetg(3,t_COMPLEX);
68 350641565 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
69 350957193 : 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 2352502 : divRc(GEN x, GEN y) {
80 2352502 : GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
81 2352494 : GEN z = cgetg(3,t_COMPLEX);
82 2352497 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
83 2352482 : gel(z,2) = gmul(mt, gel(y,2));
84 2352494 : return z;
85 : }
86 : static GEN
87 23034798 : divcR(GEN x, GEN y) {
88 23034798 : GEN z = cgetg(3,t_COMPLEX);
89 23034726 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
90 23034355 : gel(z,2) = gdiv(gel(x,2), y); return z;
91 : }
92 : static GEN
93 1274 : addRq(GEN x, GEN y) {
94 1274 : GEN z = cgetg(4,t_QUAD);
95 1274 : gel(z,1) = ZX_copy(gel(y,1));
96 1274 : gel(z,2) = gadd(x, gel(y,2));
97 1274 : gel(z,3) = gcopy(gel(y,3)); return z;
98 : }
99 : static GEN
100 2352 : mulRq(GEN x, GEN y) {
101 2352 : GEN z = cgetg(4,t_QUAD);
102 2352 : gel(z,1) = ZX_copy(gel(y,1));
103 2352 : gel(z,2) = gmul(x,gel(y,2));
104 2352 : gel(z,3) = gmul(x,gel(y,3)); return z;
105 : }
106 : static GEN
107 77 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
108 77 : long i = gexpo(x) - gexpo(y);
109 77 : if (i > 0) prec += nbits2extraprec(i);
110 77 : return gerepileupto(av, gadd(y, quadtofp(x, prec)));
111 : }
112 : static GEN
113 37415795 : mulrfrac(GEN x, GEN y)
114 : {
115 : pari_sp av;
116 37415795 : GEN z, a = gel(y,1), b = gel(y,2);
117 37415795 : if (is_pm1(a)) /* frequent special case */
118 : {
119 12618071 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
120 12619919 : return z;
121 : }
122 24797634 : av = avma;
123 24797634 : 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 4981691 : mulTp(GEN x, GEN y) { pari_sp av = avma;
151 4981691 : return gerepileupto(av, mulpp(cvtop2(x,y), y));
152 : }
153 : /* y PADIC, non zero x / y by converting x to padic */
154 : static GEN
155 3660 : divTp(GEN x, GEN y) { pari_sp av = avma;
156 3660 : return gerepileupto(av, divpp(cvtop2(x,y), y));
157 : }
158 : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
159 : * converted to O(p^e) and division by 0 */
160 : static GEN
161 1203979 : divpT(GEN x, GEN y) { pari_sp av = avma;
162 1203979 : 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 3690730 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
169 3690730 : if (lgefint(X) == 3) {
170 3367563 : ulong u = Fl_add(itou(x),itou(y), X[2]);
171 3367563 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
172 : }
173 : else {
174 323167 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
175 323166 : gel(z,2) = gerepileuptoint((pari_sp)z, u);
176 : }
177 3690730 : gel(z,1) = icopy(X); return z;
178 : }
179 : static GEN
180 1158553 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 1158553 : if (lgefint(X) == 3) {
182 784519 : ulong u = Fl_sub(itou(x),itou(y), X[2]);
183 784519 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
184 : }
185 : else {
186 374034 : GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
187 374034 : gel(z,2) = gerepileuptoint((pari_sp)z, u);
188 : }
189 1158553 : gel(z,1) = icopy(X); return z;
190 : }
191 : /* cf add_intmod_same */
192 : static GEN
193 3103529 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
194 3103529 : if (lgefint(X) == 3) {
195 2961884 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
196 2961884 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
197 : }
198 : else
199 141645 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
200 3103530 : gel(z,1) = icopy(X); return z;
201 : }
202 : /* cf add_intmod_same */
203 : static GEN
204 341865 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
205 : {
206 341865 : if (lgefint(X) == 3) {
207 319324 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
208 319317 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 22541 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
212 341858 : 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 9857930 : rfrac_denom_mul_scal(GEN d, GEN y)
225 : {
226 9857930 : GEN D = RgX_Rg_mul(d, y);
227 9857930 : 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 9857930 : return D;
233 : }
234 :
235 : static int
236 57830170 : Leading_is_neg(GEN x)
237 : {
238 122280579 : while (typ(x) == t_POL) x = leading_coeff(x);
239 57830170 : return is_real_t(typ(x))? gsigne(x) < 0: 0;
240 : }
241 :
242 : static int
243 154483885 : transtype(GEN x) { return x != gen_1 && typ(x) != t_PADIC; }
244 :
245 : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
246 : GEN
247 57834420 : gred_rfrac_simple(GEN n, GEN d)
248 : {
249 : GEN _1n, _1d, c, cn, cd, z;
250 57834420 : long dd = degpol(d);
251 :
252 57834420 : if (dd <= 0)
253 : {
254 4250 : if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
255 4250 : n = gdiv(n, gel(d,2));
256 4250 : if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
257 4250 : return n;
258 : }
259 57830170 : if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
260 57830170 : _1n = Rg_get_1(n);
261 57830170 : _1d = Rg_get_1(d);
262 57830170 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
263 57830170 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
264 57830163 : cd = content(d);
265 59700624 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
266 57830163 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
267 57830163 : if (!gequal1(cd)) {
268 6568443 : d = RgX_Rg_div(d,cd);
269 6568443 : if (!gequal1(cn))
270 : {
271 1294469 : if (gequal0(cn)) {
272 49 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
273 0 : n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
274 0 : c = gen_1;
275 : } else {
276 1294420 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
277 1294420 : c = gdiv(cn,cd);
278 : }
279 : }
280 : else
281 5273974 : c = ginv(cd);
282 : } else {
283 51261720 : if (!gequal1(cn))
284 : {
285 3277830 : if (gequal0(cn)) {
286 1484 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
287 21 : c = gen_1;
288 : } else {
289 3276346 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
290 3276346 : c = cn;
291 : }
292 : } else {
293 47983890 : GEN y = cgetg(3,t_RFRAC);
294 47983890 : gel(y,1) = gcopy(n);
295 47983890 : gel(y,2) = RgX_copy(d); return y;
296 : }
297 : }
298 :
299 9844761 : if (typ(c) == t_POL)
300 : {
301 900415 : z = c;
302 938775 : do { z = content(z); } while (typ(z) == t_POL);
303 900415 : cd = denom_i(z);
304 900415 : cn = gmul(c, cd);
305 : }
306 : else
307 : {
308 8944346 : cn = numer_i(c);
309 8944346 : cd = denom_i(c);
310 : }
311 9844761 : z = cgetg(3,t_RFRAC);
312 9844761 : gel(z,1) = gmul(n, cn);
313 9844761 : gel(z,2) = d = rfrac_denom_mul_scal(d, cd);
314 : /* can happen: Pol(O(17^-1)) / Pol([Mod(9,23), O(23^-3)]) */
315 9844761 : if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
316 9844761 : return z;
317 : }
318 :
319 : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
320 : static GEN
321 248374 : fix_rfrac(GEN x, long d)
322 : {
323 : GEN z, N, D;
324 248374 : if (!d || typ(x) == t_POL) return x;
325 171710 : z = cgetg(3, t_RFRAC);
326 171710 : N = gel(x,1);
327 171710 : D = gel(x,2);
328 171710 : if (d > 0) {
329 8722 : gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
330 165935 : : monomialcopy(N,d,varn(D));
331 165872 : gel(z, 2) = RgX_copy(D);
332 : } else {
333 5838 : gel(z, 1) = gcopy(N);
334 5838 : gel(z, 2) = RgX_shift(D, -d);
335 : }
336 171710 : return z;
337 : }
338 :
339 : /* assume d != 0 */
340 : static GEN
341 44766640 : gred_rfrac2(GEN n, GEN d)
342 : {
343 : GEN y, z, _1n, _1d;
344 : long v, vd, vn;
345 :
346 44766640 : n = simplify_shallow(n);
347 44766640 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
348 37937979 : d = simplify_shallow(d);
349 37937979 : if (typ(d) != t_POL) return gdiv(n,d);
350 36759305 : vd = varn(d);
351 36759305 : if (typ(n) != t_POL)
352 : {
353 20624450 : if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
354 20623043 : if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
355 0 : pari_err_BUG("gred_rfrac2 [incompatible variables]");
356 : }
357 16134855 : vn = varn(n);
358 16134855 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
359 15997195 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
360 15837443 : _1n = Rg_get_1(n);
361 15837443 : _1d = Rg_get_1(d);
362 15837443 : if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
363 15837443 : if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
364 :
365 : /* now n and d are t_POLs in the same variable */
366 15837443 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
367 15837443 : if (!degpol(d))
368 : {
369 12757622 : n = RgX_Rg_div(n,gel(d,2));
370 12757622 : return v? RgX_mulXn(n,v): n;
371 : }
372 :
373 : /* X does not divide gcd(n,d), deg(d) > 0 */
374 3079821 : if (!isinexact(n) && !isinexact(d))
375 : {
376 3079576 : y = RgX_divrem(n, d, &z);
377 3079576 : if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
378 248129 : z = RgX_gcd(d, z);
379 248129 : if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
380 : }
381 248374 : return fix_rfrac(gred_rfrac_simple(n,d), v);
382 : }
383 :
384 : /* x,y t_INT, return x/y in reduced form */
385 : GEN
386 134773188 : Qdivii(GEN x, GEN y)
387 : {
388 134773188 : pari_sp av = avma;
389 : GEN r, q;
390 :
391 134773188 : if (lgefint(y) == 3)
392 : {
393 116927624 : q = Qdiviu(x, y[2]);
394 116925070 : if (signe(y) > 0) return q;
395 11012100 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
396 11012267 : return q;
397 : }
398 17845564 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
399 17846984 : if (equali1(x))
400 : {
401 5178186 : if (!signe(y)) pari_err_INV("gdiv",y);
402 5178081 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
403 : }
404 12668831 : q = dvmdii(x,y,&r);
405 12668748 : if (r == gen_0) return q; /* gen_0 intended */
406 8316150 : r = gcdii(y, r);
407 8316106 : if (lgefint(r) == 3)
408 : {
409 7443053 : ulong t = r[2];
410 7443053 : set_avma(av);
411 7443014 : if (t == 1) q = mkfraccopy(x,y);
412 : else
413 : {
414 2648123 : q = cgetg(3,t_FRAC);
415 2648257 : gel(q,1) = diviuexact(x,t);
416 2648195 : gel(q,2) = diviuexact(y,t);
417 : }
418 : }
419 : else
420 : { /* rare: r and q left on stack for efficiency */
421 873053 : q = cgetg(3,t_FRAC);
422 873057 : gel(q,1) = diviiexact(x,r);
423 873057 : gel(q,2) = diviiexact(y,r);
424 : }
425 8316027 : normalize_frac(q); return q;
426 : }
427 :
428 : /* x t_INT, return x/y in reduced form */
429 : GEN
430 136237697 : Qdiviu(GEN x, ulong y)
431 : {
432 136237697 : pari_sp av = avma;
433 : ulong r, t;
434 : GEN q;
435 :
436 136237697 : if (y == 1) return icopy(x);
437 115451383 : if (!y) pari_err_INV("gdiv",gen_0);
438 115451383 : if (equali1(x)) retmkfrac(gen_1, utoipos(y));
439 108289788 : q = absdiviu_rem(x,y,&r);
440 108285665 : if (!r)
441 : {
442 62461156 : if (signe(x) < 0) togglesign(q);
443 62461237 : return q;
444 : }
445 45824509 : t = ugcd(y, r); set_avma(av);
446 45831462 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
447 18368376 : retmkfrac(diviuexact(x,t), utoipos(y / t));
448 : }
449 :
450 : /* x t_INT, return x/y in reduced form */
451 : GEN
452 1561292 : Qdivis(GEN x, long y)
453 : {
454 1561292 : pari_sp av = avma;
455 : ulong r, t;
456 : long s;
457 : GEN q;
458 :
459 1561292 : if (y > 0) return Qdiviu(x, y);
460 163500 : if (!y) pari_err_INV("gdiv",gen_0);
461 163500 : s = signe(x);
462 163500 : if (!s) return gen_0;
463 132188 : if (y < 0) { y = -y; s = -s; }
464 132188 : if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
465 131908 : if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
466 129745 : q = absdiviu_rem(x,y,&r);
467 129745 : if (!r)
468 : {
469 55356 : if (s < 0) togglesign(q);
470 55356 : return q;
471 : }
472 74389 : t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
473 74389 : if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
474 74389 : gel(q,1) = x; setsigne(x, s);
475 74389 : gel(q,2) = utoipos(y); return q;
476 : }
477 :
478 : /*******************************************************************/
479 : /* */
480 : /* CONJUGATION */
481 : /* */
482 : /*******************************************************************/
483 : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
484 : static GEN
485 18326 : quad_polmod_conj(GEN x, GEN y)
486 : {
487 : GEN z, u, v, a, b;
488 18326 : if (typ(x) != t_POL) return gcopy(x);
489 18326 : if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
490 18326 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
491 18326 : b = gel(y,3); v = gel(x,2);
492 18326 : z = cgetg(4, t_POL); z[1] = x[1];
493 18326 : gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
494 18326 : gel(z,3) = gneg(u); return z;
495 : }
496 : static GEN
497 18326 : quad_polmod_norm(GEN x, GEN y)
498 : {
499 : GEN z, u, v, a, b, c;
500 18326 : if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
501 18326 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
502 18326 : b = gel(y,3); v = gel(x,2);
503 18326 : c = gel(y,2);
504 18326 : z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
505 18326 : if (!gequal1(a)) z = gdiv(z, a);
506 18326 : return gadd(z, gsqr(v));
507 : }
508 :
509 : GEN
510 19333468 : conj_i(GEN x)
511 : {
512 : long lx, i;
513 : GEN y;
514 :
515 19333468 : switch(typ(x))
516 : {
517 6414984 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
518 6414984 : return x;
519 :
520 12755318 : case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
521 :
522 924 : case t_QUAD:
523 924 : y = cgetg(4,t_QUAD);
524 924 : gel(y,1) = gel(x,1);
525 924 : gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
526 924 : : gadd(gel(x,2), gel(x,3));
527 924 : gel(y,3) = gneg(gel(x,3)); return y;
528 :
529 9926 : case t_POL: case t_SER:
530 9926 : y = cgetg_copy(x, &lx); y[1] = x[1];
531 31913 : for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
532 9926 : return y;
533 :
534 152340 : case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
535 152340 : y = cgetg_copy(x, &lx);
536 549146 : for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
537 152338 : return y;
538 :
539 0 : case t_POLMOD:
540 : {
541 0 : GEN X = gel(x,1);
542 0 : long d = degpol(X);
543 0 : if (d < 2) return x;
544 0 : if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
545 : }
546 : }
547 0 : pari_err_TYPE("gconj",x);
548 : return NULL; /* LCOV_EXCL_LINE */
549 : }
550 : GEN
551 132418 : gconj(GEN x)
552 132418 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
553 :
554 : GEN
555 84 : conjvec(GEN x,long prec)
556 : {
557 : long lx, s, i;
558 : GEN z;
559 :
560 84 : switch(typ(x))
561 : {
562 0 : case t_INT: case t_INTMOD: case t_FRAC:
563 0 : return mkcolcopy(x);
564 :
565 0 : case t_COMPLEX: case t_QUAD:
566 0 : z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
567 :
568 28 : case t_FFELT:
569 28 : return FF_conjvec(x);
570 :
571 0 : case t_VEC: case t_COL:
572 0 : lx = lg(x); z = cgetg(lx,t_MAT);
573 0 : if (lx == 1) return z;
574 0 : gel(z,1) = conjvec(gel(x,1),prec);
575 0 : s = lgcols(z);
576 0 : for (i=2; i<lx; i++)
577 : {
578 0 : gel(z,i) = conjvec(gel(x,i),prec);
579 0 : if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
580 : }
581 0 : break;
582 :
583 56 : case t_POLMOD: {
584 56 : GEN T = gel(x,1), r;
585 : pari_sp av;
586 :
587 56 : lx = lg(T);
588 56 : if (lx <= 3) return cgetg(1,t_COL);
589 56 : x = gel(x,2);
590 238 : for (i=2; i<lx; i++)
591 : {
592 189 : GEN c = gel(T,i);
593 189 : switch(typ(c)) {
594 7 : case t_INTMOD: {
595 7 : GEN p = gel(c,1);
596 : pari_sp av;
597 7 : if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
598 7 : av = avma;
599 7 : T = RgX_to_FpX(T,p);
600 7 : x = RgX_to_FpX(x, p);
601 7 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
602 7 : z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
603 7 : return gerepileupto(av, z);
604 : }
605 182 : case t_INT:
606 182 : case t_FRAC: break;
607 0 : default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
608 : }
609 : }
610 49 : if (typ(x) != t_POL)
611 : {
612 0 : if (!is_rational_t(typ(x)))
613 0 : pari_err_TYPE("conjvec [not a rational t_POL]",x);
614 0 : retconst_col(lx-3, gcopy(x));
615 : }
616 49 : RgX_check_QX(x,"conjvec");
617 49 : av = avma;
618 49 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
619 49 : r = cleanroots(T,prec);
620 49 : z = cgetg(lx-2,t_COL);
621 182 : for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
622 49 : return gerepileupto(av, z);
623 : }
624 :
625 0 : default:
626 0 : pari_err_TYPE("conjvec",x);
627 : return NULL; /* LCOV_EXCL_LINE */
628 : }
629 0 : return z;
630 : }
631 :
632 : /********************************************************************/
633 : /** **/
634 : /** ADDITION **/
635 : /** **/
636 : /********************************************************************/
637 : /* x, y compatible PADIC, op = add or sub */
638 : static GEN
639 19953616 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
640 : {
641 19953616 : pari_sp av = avma;
642 : long d,e,r,rx,ry;
643 19953616 : GEN u, z, mod, p = gel(x,2);
644 : int swap;
645 :
646 19953616 : (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
647 19953164 : e = valp(x);
648 19953164 : r = valp(y); d = r-e;
649 19953164 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
650 19953164 : rx = precp(x);
651 19953164 : ry = precp(y);
652 19953164 : if (d) /* v(x) < v(y) */
653 : {
654 10689198 : r = d+ry; z = powiu(p,d);
655 10689333 : if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
656 10689353 : z = mulii(z,gel(y,4));
657 10689288 : u = swap? op(z, gel(x,4)): op(gel(x,4), z);
658 : }
659 : else
660 : {
661 : long c;
662 9263966 : if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
663 9263966 : u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
664 9264714 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
665 : {
666 138259 : set_avma(av); return zeropadic(p, e+r);
667 : }
668 9126316 : if (c)
669 : {
670 3351390 : mod = diviiexact(mod, powiu(p,c));
671 3351387 : r -= c;
672 3351387 : e += c;
673 : }
674 : }
675 19815584 : u = modii(u, mod);
676 19814583 : set_avma(av); z = cgetg(5,t_PADIC);
677 19813944 : z[1] = evalprecp(r) | evalvalp(e);
678 19813763 : gel(z,2) = icopy(p);
679 19813517 : gel(z,3) = icopy(mod);
680 19813269 : gel(z,4) = icopy(u); return z;
681 : }
682 : /* Rg_to_Fp(t_FRAC) without GC */
683 : static GEN
684 226838 : Q_to_Fp(GEN x, GEN p)
685 226838 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
686 : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
687 : static GEN
688 4760488 : addQp(GEN x, GEN y)
689 : {
690 4760488 : pari_sp av = avma;
691 4760488 : long d, r, e, vy = valp(y), py = precp(y);
692 4760488 : GEN z, q, mod, u, p = gel(y,2);
693 :
694 4760488 : e = Q_pvalrem(x, p, &x);
695 4760496 : d = vy - e; r = d + py;
696 4760496 : if (r <= 0) { set_avma(av); return gcopy(y); }
697 4758536 : mod = gel(y,3);
698 4758536 : u = gel(y,4);
699 4758536 : (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
700 :
701 4758526 : if (d > 0)
702 : {
703 1375930 : q = powiu(p,d);
704 1375966 : mod = mulii(mod, q);
705 1375960 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
706 1375960 : u = addii(x, mulii(u, q));
707 : }
708 3382596 : else if (d < 0)
709 : {
710 405227 : q = powiu(p,-d);
711 405227 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
712 405227 : u = addii(u, mulii(x, q));
713 405227 : r = py; e = vy;
714 : }
715 : else
716 : {
717 : long c;
718 2977369 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
719 2977369 : u = addii(u, x);
720 2977356 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
721 : {
722 1206 : set_avma(av); return zeropadic(p,e+r);
723 : }
724 2976140 : if (c)
725 : {
726 970539 : mod = diviiexact(mod, powiu(p,c));
727 970539 : r -= c;
728 970539 : e += c;
729 : }
730 : }
731 4757332 : u = modii(u, mod); set_avma(av);
732 4757295 : z = cgetg(5,t_PADIC);
733 4757246 : z[1] = evalprecp(r) | evalvalp(e);
734 4757228 : gel(z,2) = icopy(p);
735 4757238 : gel(z,3) = icopy(mod);
736 4757226 : gel(z,4) = icopy(u); return z;
737 : }
738 :
739 : /* Mod(x,X) + Mod(y,X) */
740 : #define addsub_polmod_same addsub_polmod_scal
741 : /* Mod(x,X) +/- Mod(y,Y) */
742 : static GEN
743 7203 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
744 : {
745 7203 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
746 7203 : GEN z = cgetg(3,t_POLMOD);
747 7203 : long vx = varn(X), vy = varn(Y);
748 7203 : if (vx==vy) {
749 : pari_sp av;
750 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
751 14 : warn_coercion(X,Y,gel(z,1));
752 14 : gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
753 : }
754 7189 : if (varncmp(vx, vy) < 0)
755 7189 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
756 : else
757 0 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
758 7189 : gel(z,2) = op(x, y); return z;
759 : }
760 : /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
761 : static GEN
762 13424886 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
763 13424886 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
764 :
765 : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
766 : static GEN
767 408038 : add_ser_scal(GEN y, GEN x)
768 : {
769 : long i, v, ly, vy;
770 : GEN z;
771 :
772 408038 : if (isrationalzero(x)) return gcopy(y);
773 381795 : ly = lg(y);
774 381795 : v = valser(y);
775 381795 : if (v < 3-ly) return gcopy(y);
776 : /* v + ly >= 3 */
777 381543 : if (v < 0)
778 : {
779 1162 : z = cgetg(ly,t_SER); z[1] = y[1];
780 3276 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
781 1162 : gel(z,i) = gadd(x,gel(y,i)); i++;
782 3157 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
783 1162 : return normalizeser(z);
784 : }
785 380381 : vy = varn(y);
786 380381 : if (v > 0)
787 : {
788 19271 : if (ser_isexactzero(y))
789 9310 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
790 9961 : y -= v; ly += v;
791 9961 : z = cgetg(ly,t_SER);
792 9961 : x = gcopy(x);
793 20559 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
794 : }
795 361110 : else if (ser_isexactzero(y)) return gcopy(y);
796 : else
797 : { /* v = 0, ly >= 3 */
798 361103 : z = cgetg(ly,t_SER);
799 361103 : x = gadd(x, gel(y,2));
800 361103 : i = 3;
801 : }
802 1598965 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
803 371064 : gel(z,2) = x;
804 371064 : z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
805 371064 : return gequal0(x)? normalizeser(z): z;
806 : }
807 : static long
808 7161978 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
809 : /* x,y t_SER in the same variable: x+y */
810 : static GEN
811 3581388 : ser_add(GEN x, GEN y)
812 : {
813 3581388 : long i, lx,ly, n = valser(y) - valser(x);
814 : GEN z;
815 3581388 : if (n < 0) { n = -n; swap(x,y); }
816 : /* valser(x) <= valser(y) */
817 3581388 : lx = _serprec(x);
818 3581388 : if (lx == 2) /* don't lose type information */
819 : {
820 798 : z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
821 798 : setvalser(z, valser(x)); return z;
822 : }
823 3580590 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
824 3580590 : if (n)
825 : {
826 107782 : if (n+2 > lx) return gcopy(x);
827 107082 : z = cgetg(ly,t_SER);
828 819116 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
829 506624 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
830 : } else {
831 3472808 : z = cgetg(ly,t_SER);
832 17589354 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
833 : }
834 3579890 : z[1] = x[1]; return normalizeser(z);
835 : }
836 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
837 : static GEN
838 8831008 : add_rfrac_scal(GEN y, GEN x)
839 : {
840 : pari_sp av;
841 : GEN n;
842 :
843 8831008 : if (isintzero(x)) return gcopy(y); /* frequent special case */
844 5157187 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
845 5157187 : return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
846 : }
847 :
848 : /* x "scalar", ty != t_MAT and nonscalar */
849 : static GEN
850 26104443 : add_scal(GEN y, GEN x, long ty)
851 : {
852 26104443 : switch(ty)
853 : {
854 21192447 : case t_POL: return RgX_Rg_add(y, x);
855 408010 : case t_SER: return add_ser_scal(y, x);
856 4468827 : case t_RFRAC: return add_rfrac_scal(y, x);
857 0 : case t_COL: return RgC_Rg_add(y, x);
858 35154 : case t_VEC:
859 35154 : if (isintzero(x)) return gcopy(y);
860 182 : break;
861 : }
862 187 : pari_err_TYPE2("+",x,y);
863 : return NULL; /* LCOV_EXCL_LINE */
864 : }
865 :
866 : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
867 : static GEN
868 14921043 : setfrac(GEN z, GEN a, GEN b)
869 : {
870 14921043 : gel(z,1) = icopy_avma(a, (pari_sp)z);
871 14921045 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
872 14921145 : set_avma((pari_sp)gel(z,2)); return z;
873 : }
874 : /* z <- a / (b*Q), (Q,a) = 1 */
875 : static GEN
876 14114995 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
877 : {
878 14114995 : GEN q = Qdivii(a, b); /* != 0 */
879 14115096 : if (typ(q) == t_INT)
880 : {
881 1976769 : gel(z,1) = gerepileuptoint((pari_sp)Q, q);
882 1976770 : gel(z,2) = Q; return z;
883 : }
884 12138327 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
885 : }
886 : static GEN
887 28598506 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
888 : {
889 28598506 : GEN x1 = gel(x,1), x2 = gel(x,2);
890 28598506 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
891 28598506 : int s = cmpii(x2, y2);
892 :
893 : /* common denominator: (x1 op y1) / x2 */
894 28598510 : if (!s)
895 : {
896 10144122 : pari_sp av = avma;
897 10144122 : return gerepileupto(av, Qdivii(op(x1, y1), x2));
898 : }
899 18454388 : z = cgetg(3, t_FRAC);
900 18454403 : if (s < 0)
901 : {
902 10460371 : Q = dvmdii(y2, x2, &r);
903 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
904 10460334 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
905 2762605 : delta = gcdii(x2,r);
906 : }
907 : else
908 : {
909 7994032 : Q = dvmdii(x2, y2, &r);
910 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
911 7994049 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
912 1576676 : delta = gcdii(y2,r);
913 : }
914 : /* delta = gcd(x2,y2) */
915 4339320 : if (equali1(delta))
916 : { /* numerator is nonzero */
917 1556502 : gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
918 1556501 : gel(z,2) = mulii(x2,y2); return z;
919 : }
920 2782815 : x2 = diviiexact(x2,delta);
921 2782818 : y2 = diviiexact(y2,delta);
922 2782818 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
923 2782817 : q = dvmdii(n, delta, &r);
924 2782812 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
925 2528055 : r = gcdii(delta, r);
926 2528059 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
927 2528059 : return setfrac(z, n, mulii(mulii(x2, y2), delta));
928 : }
929 :
930 : /* assume x2, y2 are t_POLs in the same variable */
931 : static GEN
932 3039726 : add_rfrac(GEN x, GEN y)
933 : {
934 3039726 : pari_sp av = avma;
935 3039726 : GEN x1 = gel(x,1), x2 = gel(x,2);
936 3039726 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
937 :
938 3039726 : delta = RgX_gcd(x2,y2);
939 3039726 : if (!degpol(delta))
940 : {
941 672 : n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
942 672 : d = RgX_mul(x2, y2);
943 672 : return gerepileupto(av, gred_rfrac_simple(n, d));
944 : }
945 3039054 : x2 = RgX_div(x2,delta);
946 3039054 : y2 = RgX_div(y2,delta);
947 3039054 : n = gadd(gmul(x1,y2), gmul(y1,x2));
948 3039054 : if (!signe(n))
949 : {
950 721425 : n = simplify_shallow(n);
951 721425 : if (isexactzero(n))
952 : {
953 721418 : if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
954 14 : return gerepileupto(av, scalarpol(n, varn(x2)));
955 : }
956 7 : return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
957 : }
958 2317629 : if (degpol(n) == 0)
959 1150596 : return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
960 1167033 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
961 1167033 : if (isexactzero(r))
962 : {
963 : GEN z;
964 228087 : d = RgX_mul(x2, y2);
965 : /* "constant" denominator ? */
966 228087 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
967 228087 : return gerepileupto(av, z);
968 : }
969 938946 : r = RgX_gcd(delta, r);
970 938946 : if (degpol(r))
971 : {
972 160740 : n = RgX_div(n, r);
973 160740 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
974 : }
975 : else
976 778206 : d = RgX_mul(gel(x,2), y2);
977 938946 : return gerepileupto(av, gred_rfrac_simple(n, d));
978 : }
979 :
980 : GEN
981 5882638098 : gadd(GEN x, GEN y)
982 : {
983 5882638098 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
984 : pari_sp av, tetpil;
985 : GEN z, p1;
986 :
987 5882638098 : if (tx == ty) switch(tx) /* shortcut to generic case */
988 : {
989 2911263149 : case t_INT: return addii(x,y);
990 1636380362 : case t_REAL: return addrr(x,y);
991 1521962 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
992 1521962 : z = cgetg(3,t_INTMOD);
993 1521962 : if (X==Y || equalii(X,Y))
994 1521948 : return add_intmod_same(z, X, gel(x,2), gel(y,2));
995 14 : gel(z,1) = gcdii(X,Y);
996 14 : warn_coercion(X,Y,gel(z,1));
997 14 : av = avma; p1 = addii(gel(x,2),gel(y,2));
998 14 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
999 : }
1000 23164566 : case t_FRAC: return addsub_frac(x,y,addii);
1001 354178174 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1002 353749592 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1003 354550120 : if (isintzero(gel(z,2)))
1004 : {
1005 521178 : set_avma((pari_sp)(z+3));
1006 521178 : return gadd(gel(x,1),gel(y,1));
1007 : }
1008 353950901 : gel(z,1) = gadd(gel(x,1),gel(y,1));
1009 353961315 : return z;
1010 14516033 : case t_PADIC:
1011 14516033 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
1012 14515822 : return addsub_pp(x,y, addii);
1013 672 : case t_QUAD: z = cgetg(4,t_QUAD);
1014 672 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1015 672 : gel(z,1) = ZX_copy(gel(x,1));
1016 672 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1017 672 : gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
1018 10161036 : case t_POLMOD:
1019 10161036 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1020 10153868 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
1021 7168 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
1022 14006799 : case t_FFELT: return FF_add(x,y);
1023 72883985 : case t_POL:
1024 72883985 : vx = varn(x);
1025 72883985 : vy = varn(y);
1026 72883985 : if (vx != vy) {
1027 852726 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1028 65287 : else return RgX_Rg_add(y, x);
1029 : }
1030 72031259 : return RgX_add(x, y);
1031 3576033 : case t_SER:
1032 3576033 : vx = varn(x);
1033 3576033 : vy = varn(y);
1034 3576033 : if (vx != vy) {
1035 28 : if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1036 21 : else return add_ser_scal(y, x);
1037 : }
1038 3576005 : return ser_add(x, y);
1039 4331627 : case t_RFRAC:
1040 4331627 : vx = varn(gel(x,2));
1041 4331627 : vy = varn(gel(y,2));
1042 4331627 : if (vx != vy) {
1043 1291901 : if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1044 538397 : else return add_rfrac_scal(y, x);
1045 : }
1046 3039726 : return add_rfrac(x,y);
1047 2035466 : case t_VEC:
1048 2035466 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1049 2035466 : return RgV_add(x,y);
1050 1101620 : case t_COL:
1051 1101620 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1052 1101620 : return RgC_add(x,y);
1053 671356 : case t_MAT:
1054 671356 : lx = lg(x);
1055 671356 : if (lg(y) != lx) pari_err_OP("+",x,y);
1056 671356 : if (lx == 1) return cgetg(1, t_MAT);
1057 671356 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1058 671349 : return RgM_add(x,y);
1059 :
1060 0 : default: pari_err_TYPE2("+",x,y);
1061 : }
1062 : /* tx != ty */
1063 837130882 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1064 :
1065 837130882 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1066 : {
1067 740430035 : case t_INT:
1068 : switch(ty)
1069 : {
1070 436510253 : case t_REAL: return addir(x,y);
1071 2139209 : case t_INTMOD:
1072 2139209 : z = cgetg(3, t_INTMOD);
1073 2139209 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1074 45009667 : case t_FRAC: z = cgetg(3,t_FRAC);
1075 45009664 : gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1076 45009618 : gel(z,2) = icopy(gel(y,2)); return z;
1077 249859751 : case t_COMPLEX: return addRc(x, y);
1078 5594400 : case t_PADIC:
1079 5594400 : if (!signe(x)) return gcopy(y);
1080 4530752 : return addQp(x,y);
1081 1113 : case t_QUAD: return addRq(x, y);
1082 1357786 : case t_FFELT: return FF_Z_add(y,x);
1083 : }
1084 :
1085 : case t_REAL:
1086 57076523 : switch(ty)
1087 : {
1088 13518545 : case t_FRAC:
1089 13518545 : if (!signe(gel(y,1))) return rcopy(x);
1090 13518545 : if (!signe(x))
1091 : {
1092 11865 : lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1093 11865 : return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1094 : }
1095 13506680 : av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
1096 13505079 : return gerepile(av,tetpil,divri(z,gel(y,2)));
1097 43557908 : case t_COMPLEX: return addRc(x, y);
1098 70 : case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, realprec(x));
1099 :
1100 0 : default: pari_err_TYPE2("+",x,y);
1101 : }
1102 :
1103 17640 : case t_INTMOD:
1104 : switch(ty)
1105 : {
1106 17507 : case t_FRAC: { GEN X = gel(x,1);
1107 17507 : z = cgetg(3, t_INTMOD);
1108 17507 : p1 = Fp_div(gel(y,1), gel(y,2), X);
1109 17507 : return add_intmod_same(z, X, p1, gel(x,2));
1110 : }
1111 14 : case t_FFELT:
1112 14 : if (!equalii(gel(x,1),FF_p_i(y)))
1113 0 : pari_err_OP("+",x,y);
1114 14 : return FF_Z_add(y,gel(x,2));
1115 91 : case t_COMPLEX: return addRc(x, y);
1116 0 : case t_PADIC: { GEN X = gel(x,1);
1117 0 : z = cgetg(3, t_INTMOD);
1118 0 : return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1119 : }
1120 28 : case t_QUAD: return addRq(x, y);
1121 : }
1122 :
1123 : case t_FRAC:
1124 : switch (ty)
1125 : {
1126 10296597 : case t_COMPLEX: return addRc(x, y);
1127 229736 : case t_PADIC:
1128 229736 : if (!signe(gel(x,1))) return gcopy(y);
1129 229736 : return addQp(x,y);
1130 133 : case t_QUAD: return addRq(x, y);
1131 1337 : case t_FFELT: return FF_Q_add(y, x);
1132 : }
1133 :
1134 : case t_FFELT:
1135 1 : pari_err_TYPE2("+",x,y);
1136 :
1137 35 : case t_COMPLEX:
1138 : switch(ty)
1139 : {
1140 28 : case t_PADIC:
1141 28 : return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
1142 7 : case t_QUAD:
1143 7 : lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1144 7 : return gequal0(y)? gcopy(x): addqf(y, x, lx);
1145 : }
1146 :
1147 : case t_PADIC: /* ty == t_QUAD */
1148 0 : return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
1149 : }
1150 : /* tx < ty, !is_const_t(y) */
1151 31815490 : switch(ty)
1152 : {
1153 7951 : case t_MAT:
1154 7951 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1155 7951 : if (isrationalzero(x)) return gcopy(y);
1156 7874 : return RgM_Rg_add(y, x);
1157 206924 : case t_COL:
1158 206924 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1159 206924 : return RgC_Rg_add(y, x);
1160 2392180 : case t_POLMOD: /* is_const_t(tx) in this case */
1161 2392180 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1162 : }
1163 29208435 : if (is_scalar_t(tx)) {
1164 26109542 : if (tx == t_POLMOD)
1165 : {
1166 112774 : vx = varn(gel(x,1));
1167 112774 : vy = gvar(y);
1168 112773 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1169 : else
1170 85487 : if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1171 27811 : return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1172 : }
1173 25996768 : return add_scal(y, x, ty);
1174 : }
1175 : /* x and y are not scalars, ty != t_MAT */
1176 3098903 : vx = gvar(x);
1177 3098905 : vy = gvar(y);
1178 3098905 : if (vx != vy) { /* x or y is treated as a scalar */
1179 22717 : if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1180 32340 : return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1181 32340 : : add_scal(y, x, ty);
1182 : }
1183 : /* vx = vy */
1184 3076188 : switch(tx)
1185 : {
1186 3075691 : case t_POL:
1187 : switch (ty)
1188 : {
1189 5411 : case t_SER:
1190 5411 : if (lg(x) == 2) return gcopy(y);
1191 5390 : i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1192 5390 : i = lg(y) + valser(y) - i;
1193 5390 : if (i < 3) return gcopy(y);
1194 5383 : p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1195 5383 : settyp(p1, t_VECSMALL); /* p1 left on stack */
1196 5383 : return y;
1197 :
1198 3070280 : case t_RFRAC: return add_rfrac_scal(y, x);
1199 : }
1200 0 : break;
1201 :
1202 497 : case t_SER:
1203 497 : if (ty == t_RFRAC)
1204 : {
1205 : long vn, vd;
1206 497 : av = avma;
1207 497 : vn = gval(gel(y,1), vy);
1208 497 : vd = RgX_valrem_inexact(gel(y,2), NULL);
1209 497 : if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
1210 :
1211 497 : l = lg(x) + valser(x) - (vn - vd);
1212 497 : if (l < 3) { set_avma(av); return gcopy(x); }
1213 497 : return gerepileupto(av, gadd(x, rfrac_to_ser_i(y, l)));
1214 : }
1215 0 : break;
1216 : }
1217 0 : pari_err_TYPE2("+",x,y);
1218 : return NULL; /* LCOV_EXCL_LINE */
1219 : }
1220 :
1221 : GEN
1222 270554594 : gaddsg(long x, GEN y)
1223 : {
1224 270554594 : long ty = typ(y);
1225 : GEN z;
1226 :
1227 270554594 : switch(ty)
1228 : {
1229 124828776 : case t_INT: return addsi(x,y);
1230 119174876 : case t_REAL: return addsr(x,y);
1231 12011 : case t_INTMOD:
1232 12011 : z = cgetg(3, t_INTMOD);
1233 12011 : return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1234 15247397 : case t_FRAC: z = cgetg(3,t_FRAC);
1235 15247397 : gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1236 15247397 : gel(z,2) = icopy(gel(y,2)); return z;
1237 8163832 : case t_COMPLEX:
1238 8163832 : z = cgetg(3, t_COMPLEX);
1239 8163832 : gel(z,1) = gaddsg(x, gel(y,1));
1240 8163832 : gel(z,2) = gcopy(gel(y,2)); return z;
1241 :
1242 3127702 : default: return gadd(stoi(x), y);
1243 : }
1244 : }
1245 :
1246 : GEN
1247 3166761 : gsubsg(long x, GEN y)
1248 : {
1249 : GEN z, a, b;
1250 : pari_sp av;
1251 :
1252 3166761 : switch(typ(y))
1253 : {
1254 276639 : case t_INT: return subsi(x,y);
1255 1224645 : case t_REAL: return subsr(x,y);
1256 56 : case t_INTMOD:
1257 56 : z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1258 56 : return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1259 732302 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1260 732302 : gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
1261 732302 : gel(z,2) = icopy(gel(y,2)); return z;
1262 896005 : case t_COMPLEX:
1263 896005 : z = cgetg(3, t_COMPLEX);
1264 896005 : gel(z,1) = gsubsg(x, gel(y,1));
1265 896005 : gel(z,2) = gneg(gel(y,2)); return z;
1266 : }
1267 37114 : av = avma;
1268 37114 : return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
1269 : }
1270 :
1271 : /********************************************************************/
1272 : /** **/
1273 : /** SUBTRACTION **/
1274 : /** **/
1275 : /********************************************************************/
1276 :
1277 : GEN
1278 2827922550 : gsub(GEN x, GEN y)
1279 : {
1280 2827922550 : long tx = typ(x), ty = typ(y);
1281 : pari_sp av;
1282 : GEN z;
1283 2827922550 : if (tx == ty) switch(tx) /* shortcut to generic case */
1284 : {
1285 2059621860 : case t_INT: return subii(x,y);
1286 576778761 : case t_REAL: return subrr(x,y);
1287 1158567 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1288 1158567 : z = cgetg(3,t_INTMOD);
1289 1158567 : if (X==Y || equalii(X,Y))
1290 1158553 : return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1291 14 : gel(z,1) = gcdii(X,Y);
1292 14 : warn_coercion(X,Y,gel(z,1));
1293 14 : av = avma; p1 = subii(gel(x,2),gel(y,2));
1294 14 : gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
1295 : }
1296 5433949 : case t_FRAC: return addsub_frac(x,y, subii);
1297 101903938 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1298 101944199 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1299 101849046 : if (isintzero(gel(z,2)))
1300 : {
1301 21301 : set_avma((pari_sp)(z+3));
1302 21301 : return gsub(gel(x,1),gel(y,1));
1303 : }
1304 101832257 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1305 101814073 : return z;
1306 5437800 : case t_PADIC:
1307 5437800 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
1308 5437800 : return addsub_pp(x,y, subii);
1309 168 : case t_QUAD: z = cgetg(4,t_QUAD);
1310 168 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1311 168 : gel(z,1) = ZX_copy(gel(x,1));
1312 168 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1313 168 : gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1314 851062 : case t_POLMOD:
1315 851062 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1316 851027 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1317 35 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1318 203215 : case t_FFELT: return FF_sub(x,y);
1319 10082420 : case t_POL: {
1320 10082420 : long vx = varn(x);
1321 10082420 : long vy = varn(y);
1322 10082420 : if (vx != vy) {
1323 22953 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1324 1988 : else return Rg_RgX_sub(x, y);
1325 : }
1326 10059467 : return RgX_sub(x, y);
1327 : }
1328 299247 : case t_VEC:
1329 299247 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1330 299247 : return RgV_sub(x,y);
1331 3567141 : case t_COL:
1332 3567141 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1333 3567141 : return RgC_sub(x,y);
1334 251830 : case t_MAT: {
1335 251830 : long lx = lg(x);
1336 251830 : if (lg(y) != lx) pari_err_OP("+",x,y);
1337 251831 : if (lx == 1) return cgetg(1, t_MAT);
1338 172260 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1339 172260 : return RgM_sub(x,y);
1340 : }
1341 2469821 : case t_RFRAC: case t_SER: break;
1342 :
1343 0 : default: pari_err_TYPE2("+",x,y);
1344 : }
1345 63980358 : av = avma;
1346 63980358 : return gerepileupto(av, gadd(x,gneg_i(y)));
1347 : }
1348 :
1349 : /********************************************************************/
1350 : /** **/
1351 : /** MULTIPLICATION **/
1352 : /** **/
1353 : /********************************************************************/
1354 : static GEN
1355 316204 : mul_ser_scal(GEN y, GEN x)
1356 : {
1357 : long l, i;
1358 : GEN z;
1359 316204 : if (isexactzero(x)) return gmul(Rg_get_0(y), x);
1360 312956 : if (ser_isexactzero(y))
1361 : {
1362 525 : z = scalarser(lg(y) == 2? Rg_get_0(x): gmul(gel(y,2), x), varn(y), 1);
1363 525 : setvalser(z, valser(y)); return z;
1364 : }
1365 312431 : z = cgetg_copy(y, &l); z[1] = y[1];
1366 4218550 : for (i = 2; i < l; i++) gel(z,i) = gmul(gel(y,i), x);
1367 312431 : return normalizeser(z);
1368 : }
1369 : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1370 : * [n/d a valid RFRAC] */
1371 : static GEN
1372 10453385 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1373 : {
1374 10453385 : pari_sp av = avma;
1375 : GEN z;
1376 :
1377 10453385 : switch(typ(x))
1378 : {
1379 21 : case t_PADIC:
1380 21 : n = gmul(n, x);
1381 21 : d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
1382 21 : return gerepileupto(av, gdiv(n,d));
1383 :
1384 518 : case t_INTMOD: case t_POLMOD:
1385 518 : n = gmul(n, x);
1386 518 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1387 518 : return gerepileupto(av, gdiv(n,d));
1388 : }
1389 10452846 : z = gred_rfrac2(x, d);
1390 10452846 : n = simplify_shallow(n);
1391 10452846 : if (typ(z) == t_RFRAC)
1392 : {
1393 7936833 : n = gmul(gel(z,1), n);
1394 7936833 : d = gel(z,2);
1395 7936833 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1396 0 : z = RgX_Rg_div(n, d);
1397 : else
1398 7936833 : z = gred_rfrac_simple(n, d);
1399 : }
1400 : else
1401 2516013 : z = gmul(z, n);
1402 10452846 : return gerepileupto(av, z);
1403 : }
1404 : static GEN
1405 121949729 : mul_scal(GEN y, GEN x, long ty)
1406 : {
1407 121949729 : switch(ty)
1408 : {
1409 113054885 : case t_POL:
1410 113054885 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1411 109904933 : return RgX_Rg_mul(y, x);
1412 308063 : case t_SER: return mul_ser_scal(y, x);
1413 8586788 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1414 14 : case t_QFB:
1415 14 : if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1416 : }
1417 12 : pari_err_TYPE2("*",x,y);
1418 : return NULL; /* LCOV_EXCL_LINE */
1419 : }
1420 :
1421 : static GEN
1422 160865 : mul_gen_rfrac(GEN X, GEN Y)
1423 : {
1424 160865 : GEN y1 = gel(Y,1), y2 = gel(Y,2);
1425 160865 : long vx = gvar(X), vy = varn(y2);
1426 166626 : return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1427 5761 : gred_rfrac_simple(gmul(y1,X), y2);
1428 : }
1429 : /* (x1/x2) * (y1/y2) */
1430 : static GEN
1431 7907135 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1432 : {
1433 : GEN z, X, Y;
1434 7907135 : pari_sp av = avma;
1435 :
1436 7907135 : X = gred_rfrac2(x1, y2);
1437 7907135 : Y = gred_rfrac2(y1, x2);
1438 7907135 : if (typ(X) == t_RFRAC)
1439 : {
1440 6628053 : if (typ(Y) == t_RFRAC) {
1441 6562477 : x1 = gel(X,1);
1442 6562477 : x2 = gel(X,2);
1443 6562477 : y1 = gel(Y,1);
1444 6562477 : y2 = gel(Y,2);
1445 6562477 : z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1446 : } else
1447 65576 : z = mul_gen_rfrac(Y, X);
1448 : }
1449 1279082 : else if (typ(Y) == t_RFRAC)
1450 95289 : z = mul_gen_rfrac(X, Y);
1451 : else
1452 1183793 : z = gmul(X, Y);
1453 7907135 : return gerepileupto(av, z);
1454 : }
1455 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1456 : static GEN
1457 271003 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1458 : {
1459 271003 : pari_sp av = avma;
1460 271003 : GEN X = gred_rfrac2(x1, y2);
1461 271003 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1462 : {
1463 264423 : x2 = RgX_mul(gel(X,2), x2);
1464 264423 : x1 = gel(X,1);
1465 : }
1466 : else
1467 6580 : x1 = X;
1468 271003 : return gerepileupto(av, gred_rfrac_simple(x1, x2));
1469 : }
1470 :
1471 : /* Mod(y, Y) * x, assuming x scalar */
1472 : static GEN
1473 3427483 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1474 3427483 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1475 :
1476 : /* cf mulqq */
1477 : static GEN
1478 5858586 : quad_polmod_mul(GEN P, GEN x, GEN y)
1479 : {
1480 5858586 : GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
1481 5858586 : pari_sp tetpil, av = avma;
1482 5858586 : T[1] = x[1];
1483 5858586 : p2 = gmul(gel(x,2), gel(y,2));
1484 5858586 : p3 = gmul(gel(x,3), gel(y,3));
1485 5858586 : p1 = gmul(gneg_i(c),p3);
1486 : /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
1487 5858586 : if (typ(b) == t_INT)
1488 : {
1489 5858565 : if (signe(b))
1490 : {
1491 4284101 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1492 4284101 : if (is_pm1(b))
1493 : {
1494 4283443 : if (signe(b) > 0) p3 = gneg(p3);
1495 : }
1496 : else
1497 658 : p3 = gmul(negi(b), p3);
1498 : }
1499 : else
1500 : {
1501 1574464 : p3 = gmul(gel(x,2),gel(y,3));
1502 1574464 : p4 = gmul(gel(x,3),gel(y,2));
1503 : }
1504 : }
1505 : else
1506 : {
1507 21 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1508 21 : p3 = gmul(gneg_i(b), p3);
1509 : }
1510 5858586 : tetpil = avma;
1511 5858586 : gel(T,2) = gadd(p2, p1);
1512 5858586 : gel(T,3) = gadd(p4, p3);
1513 5858586 : gerepilecoeffssp(av,tetpil,T+2,2);
1514 5858586 : return normalizepol_lg(T,4);
1515 : }
1516 : /* Mod(x,T) * Mod(y,T) */
1517 : static GEN
1518 8562223 : mul_polmod_same(GEN T, GEN x, GEN y)
1519 : {
1520 8562223 : GEN z = cgetg(3,t_POLMOD), a;
1521 8562223 : long v = varn(T), lx = lg(x), ly = lg(y);
1522 8562223 : gel(z,1) = RgX_copy(T);
1523 : /* x * y mod T optimised */
1524 8562223 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1525 7891506 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1526 1587130 : a = gmul(x, y);
1527 : else
1528 : {
1529 6975093 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1530 5853315 : a = quad_polmod_mul(T, x, y);
1531 : else
1532 1121778 : a = RgXQ_mul(x, y, gel(z,1));
1533 : }
1534 8562223 : gel(z,2) = a; return z;
1535 : }
1536 : static GEN
1537 484324 : sqr_polmod(GEN T, GEN x)
1538 : {
1539 484324 : GEN a, z = cgetg(3,t_POLMOD);
1540 484324 : gel(z,1) = RgX_copy(T);
1541 484324 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1542 93668 : a = gsqr(x);
1543 : else
1544 : {
1545 390656 : pari_sp av = avma;
1546 390656 : a = RgXQ_sqr(x, gel(z,1));
1547 390656 : a = gerepileupto(av, a);
1548 : }
1549 484324 : gel(z,2) = a; return z;
1550 : }
1551 : /* Mod(x,X) * Mod(y,Y) */
1552 : static GEN
1553 2675106 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1554 : {
1555 2675106 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1556 2675106 : long vx = varn(X), vy = varn(Y);
1557 2675106 : GEN z = cgetg(3,t_POLMOD);
1558 :
1559 2675106 : if (vx==vy) {
1560 : pari_sp av;
1561 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
1562 14 : warn_coercion(X,Y,gel(z,1));
1563 14 : gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
1564 14 : return z;
1565 : }
1566 2675092 : if (varncmp(vx, vy) < 0)
1567 410662 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1568 : else
1569 2264430 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1570 2675092 : gel(z,2) = gmul(x, y); return z;
1571 : }
1572 :
1573 : #if 0 /* used by 3M only */
1574 : /* set z = x+y and return 1 if x,y have the same sign
1575 : * set z = x-y and return 0 otherwise */
1576 : static int
1577 : did_add(GEN x, GEN y, GEN *z)
1578 : {
1579 : long tx = typ(x), ty = typ(y);
1580 : if (tx == ty) switch(tx)
1581 : {
1582 : case t_INT: *z = addii(x,y); return 1;
1583 : case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1584 : case t_REAL:
1585 : if (signe(x) == -signe(y))
1586 : { *z = subrr(x,y); return 0; }
1587 : else
1588 : { *z = addrr(x,y); return 1; }
1589 : }
1590 : if (tx == t_REAL) switch(ty)
1591 : {
1592 : case t_INT:
1593 : if (signe(x) == -signe(y))
1594 : { *z = subri(x,y); return 0; }
1595 : else
1596 : { *z = addri(x,y); return 1; }
1597 : case t_FRAC:
1598 : if (signe(x) == -signe(gel(y,1)))
1599 : { *z = gsub(x,y); return 0; }
1600 : else
1601 : { *z = gadd(x,y); return 1; }
1602 : }
1603 : else if (ty == t_REAL) switch(tx)
1604 : {
1605 : case t_INT:
1606 : if (signe(x) == -signe(y))
1607 : { *z = subir(x,y); return 0; }
1608 : else
1609 : { *z = addir(x,y); return 1; }
1610 : case t_FRAC:
1611 : if (signe(gel(x,1)) == -signe(y))
1612 : { *z = gsub(x,y); return 0; }
1613 : else
1614 : { *z = gadd(x,y); return 1; }
1615 : }
1616 : *z = gadd(x,y); return 1;
1617 : }
1618 : #endif
1619 : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1620 : static GEN
1621 11645643 : mulcIR(GEN x, GEN y)
1622 : {
1623 11645643 : GEN z = cgetg(3,t_COMPLEX);
1624 11645776 : pari_sp av = avma;
1625 11645776 : gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
1626 11645876 : gel(z,2) = gmul(y, gel(x,1));
1627 11645847 : return z;
1628 :
1629 : }
1630 : /* x,y COMPLEX */
1631 : static GEN
1632 273441623 : mulcc(GEN x, GEN y)
1633 : {
1634 273441623 : GEN xr = gel(x,1), xi = gel(x,2);
1635 273441623 : GEN yr = gel(y,1), yi = gel(y,2);
1636 : GEN p1, p2, p3, p4, z;
1637 : pari_sp tetpil, av;
1638 :
1639 273441623 : if (isintzero(xr))
1640 : {
1641 15627496 : if (isintzero(yr)) {
1642 7122562 : av = avma;
1643 7122562 : return gerepileupto(av, gneg(gmul(xi,yi)));
1644 : }
1645 8504773 : return mulcIR(y, xi);
1646 : }
1647 257807061 : if (isintzero(yr)) return mulcIR(x, yi);
1648 :
1649 254668285 : z = cgetg(3,t_COMPLEX); av = avma;
1650 : #if 0
1651 : /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1652 : * e.g. xr + xi if exponents differ */
1653 : if (did_add(xr, xi, &p3))
1654 : {
1655 : if (did_add(yr, yi, &p4)) {
1656 : /* R = xr*yr - xi*yi
1657 : * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1658 : p1 = gmul(xr,yr);
1659 : p2 = gmul(xi,yi); p2 = gneg(p2);
1660 : p3 = gmul(p3, p4);
1661 : p4 = gsub(p2, p1);
1662 : } else {
1663 : /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1664 : * I = xr*yi + xi*yr */
1665 : p1 = gmul(p3,p4);
1666 : p3 = gmul(xr,yi);
1667 : p4 = gmul(xi,yr);
1668 : p2 = gsub(p3, p4);
1669 : }
1670 : } else {
1671 : if (did_add(yr, yi, &p4)) {
1672 : /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1673 : * I = xr*yi +xi*yr */
1674 : p1 = gmul(p3,p4);
1675 : p3 = gmul(xr,yi);
1676 : p4 = gmul(xi,yr);
1677 : p2 = gsub(p4, p3);
1678 : } else {
1679 : /* R = xr*yr - xi*yi
1680 : * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1681 : p3 = gneg( gmul(p3, p4) );
1682 : p1 = gmul(xr,yr);
1683 : p2 = gmul(xi,yi);
1684 : p4 = gadd(p1, p2);
1685 :
1686 : p2 = gneg(p2);
1687 : }
1688 : }
1689 : tetpil = avma;
1690 : gel(z,1) = gadd(p1,p2);
1691 : gel(z,2) = gadd(p3,p4);
1692 : #else
1693 254657653 : if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1694 : { /* 3M formula */
1695 559286 : p3 = addii(xr,xi);
1696 559286 : p4 = addii(yr,yi);
1697 559286 : p1 = mulii(xr,yr);
1698 559286 : p2 = mulii(xi,yi);
1699 559286 : p3 = mulii(p3,p4);
1700 559286 : p4 = addii(p2,p1);
1701 559286 : tetpil = avma;
1702 559286 : gel(z,1) = subii(p1,p2);
1703 559286 : gel(z,2) = subii(p3,p4);
1704 559286 : if (!signe(gel(z,2)))
1705 113225 : return gerepileuptoint((pari_sp)(z+3), gel(z,1));
1706 : }
1707 : else
1708 : { /* naive 4M formula: avoid all loss of accuracy */
1709 254098367 : p1 = gmul(xr,yr);
1710 254027698 : p2 = gmul(xi,yi);
1711 254029066 : p3 = gmul(xr,yi);
1712 254034466 : p4 = gmul(xi,yr);
1713 254042113 : tetpil = avma;
1714 254042113 : gel(z,1) = gsub(p1,p2);
1715 253912986 : gel(z,2) = gadd(p3,p4);
1716 253916585 : if (isintzero(gel(z,2)))
1717 : {
1718 50575 : cgiv(gel(z,2));
1719 50575 : return gerepileupto((pari_sp)(z+3), gel(z,1));
1720 : }
1721 : }
1722 : #endif
1723 :
1724 254302259 : gerepilecoeffssp(av,tetpil, z+1,2); return z;
1725 : }
1726 : /* x,y PADIC */
1727 : static GEN
1728 17470418 : mulpp(GEN x, GEN y) {
1729 17470418 : long l = valp(x) + valp(y);
1730 : pari_sp av;
1731 : GEN z, t;
1732 17470418 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
1733 17470391 : if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
1734 17243737 : if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
1735 :
1736 16897190 : t = (precp(x) > precp(y))? y: x;
1737 16897190 : z = cgetp(t); setvalp(z,l); av = avma;
1738 16897277 : affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
1739 16897106 : set_avma(av); return z;
1740 : }
1741 : /* x,y QUAD */
1742 : static GEN
1743 1106 : mulqq(GEN x, GEN y) {
1744 1106 : GEN z = cgetg(4,t_QUAD);
1745 1106 : GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
1746 : pari_sp av, tetpil;
1747 1106 : if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
1748 :
1749 1106 : gel(z,1) = ZX_copy(P); av = avma;
1750 1106 : p2 = gmul(gel(x,2),gel(y,2));
1751 1106 : p3 = gmul(gel(x,3),gel(y,3));
1752 1106 : p1 = gmul(gneg_i(c),p3);
1753 :
1754 1106 : if (signe(b))
1755 987 : p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
1756 : else
1757 : {
1758 119 : p3 = gmul(gel(x,2),gel(y,3));
1759 119 : p4 = gmul(gel(x,3),gel(y,2));
1760 : }
1761 1106 : tetpil = avma;
1762 1106 : gel(z,2) = gadd(p2,p1);
1763 1106 : gel(z,3) = gadd(p4,p3);
1764 1106 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
1765 : }
1766 :
1767 : GEN
1768 15544110 : mulcxI(GEN x)
1769 : {
1770 : GEN z;
1771 15544110 : switch(typ(x))
1772 : {
1773 2627988 : case t_INT: case t_REAL: case t_FRAC:
1774 2627988 : return mkcomplex(gen_0, x);
1775 12916639 : case t_COMPLEX:
1776 12916639 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1777 12903214 : z = cgetg(3,t_COMPLEX);
1778 12903576 : gel(z,1) = gneg(gel(x,2));
1779 12904045 : gel(z,2) = gel(x,1); return z;
1780 48 : default:
1781 48 : return gmul(gen_I(), x);
1782 : }
1783 : }
1784 : GEN
1785 3060948 : mulcxmI(GEN x)
1786 : {
1787 : GEN z;
1788 3060948 : switch(typ(x))
1789 : {
1790 116059 : case t_INT: case t_REAL: case t_FRAC:
1791 116059 : return mkcomplex(gen_0, gneg(x));
1792 2760163 : case t_COMPLEX:
1793 2760163 : if (isintzero(gel(x,1))) return gel(x,2);
1794 2721983 : z = cgetg(3,t_COMPLEX);
1795 2721964 : gel(z,1) = gel(x,2);
1796 2721964 : gel(z,2) = gneg(gel(x,1)); return z;
1797 184726 : default:
1798 184726 : return gmul(mkcomplex(gen_0, gen_m1), x);
1799 : }
1800 : }
1801 : /* x * I^k */
1802 : GEN
1803 5522716 : mulcxpowIs(GEN x, long k)
1804 : {
1805 5522716 : switch (k & 3)
1806 : {
1807 1368890 : case 1: return mulcxI(x);
1808 1361896 : case 2: return gneg(x);
1809 1363488 : case 3: return mulcxmI(x);
1810 : }
1811 1428442 : return x;
1812 : }
1813 :
1814 : /* caller will assume l > 2 later */
1815 : static GEN
1816 5706179 : init_ser(long l, long v, long e)
1817 : {
1818 5706179 : GEN z = cgetg(l, t_SER);
1819 5706179 : z[1] = evalvalser(e) | evalvarn(v) | evalsigne(1); return z;
1820 : }
1821 :
1822 : /* fill in coefficients of t_SER z from coeffs of t_POL y */
1823 : static GEN
1824 5693846 : fill_ser(GEN z, GEN y)
1825 : {
1826 5693846 : long i, lx = lg(z), ly = lg(y); /* lx > 2 */
1827 5693846 : if (ly >= lx) {
1828 25570417 : for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
1829 : } else {
1830 335415 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1831 239573 : for ( ; i < lx; i++) gel(z,i) = gen_0;
1832 : }
1833 : /* dangerous case (already normalized), don't use normalizeser */
1834 5693846 : if (ly == 3 && !signe(y)) { setsigne(z, 0); return z; }
1835 5693748 : return normalizeser(z);
1836 : }
1837 : /* assume typ(x) = t_VEC */
1838 : static int
1839 98 : is_ext_qfr(GEN x)
1840 91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
1841 189 : && typ(gel(x,2)) == t_REAL; }
1842 :
1843 : GEN
1844 8515630828 : gmul(GEN x, GEN y)
1845 : {
1846 : long tx, ty, lx, ly, vx, vy, i, l;
1847 : pari_sp av, tetpil;
1848 : GEN z, p1;
1849 :
1850 8515630828 : if (x == y) return gsqr(x);
1851 7637511243 : tx = typ(x); ty = typ(y);
1852 7637511243 : if (tx == ty) switch(tx)
1853 : {
1854 3623425559 : case t_INT: return mulii(x,y);
1855 1866173181 : case t_REAL: return mulrr(x,y);
1856 1780384 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1857 1780384 : z = cgetg(3,t_INTMOD);
1858 1780384 : if (X==Y || equalii(X,Y))
1859 1780349 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1860 35 : gel(z,1) = gcdii(X,Y);
1861 35 : warn_coercion(X,Y,gel(z,1));
1862 35 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1863 35 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
1864 : }
1865 23410303 : case t_FRAC:
1866 : {
1867 23410303 : GEN x1 = gel(x,1), x2 = gel(x,2);
1868 23410303 : GEN y1 = gel(y,1), y2 = gel(y,2);
1869 23410303 : z=cgetg(3,t_FRAC);
1870 23410305 : p1 = gcdii(x1, y2);
1871 23410301 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1872 23410301 : p1 = gcdii(x2, y1);
1873 23410300 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1874 23410300 : tetpil = avma;
1875 23410300 : gel(z,2) = mulii(x2,y2);
1876 23410302 : gel(z,1) = mulii(x1,y1);
1877 23410303 : fix_frac_if_int_GC(z,tetpil); return z;
1878 : }
1879 264485873 : case t_COMPLEX: return mulcc(x, y);
1880 12480097 : case t_PADIC: return mulpp(x, y);
1881 882 : case t_QUAD: return mulqq(x, y);
1882 14546045 : case t_FFELT: return FF_mul(x, y);
1883 11002430 : case t_POLMOD:
1884 11002430 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1885 8327324 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1886 2675106 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1887 48061630 : case t_POL:
1888 48061630 : vx = varn(x);
1889 48061630 : vy = varn(y);
1890 48061630 : if (vx != vy) {
1891 4847531 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1892 2067116 : else return RgX_Rg_mul(y, x);
1893 : }
1894 43214099 : return RgX_mul(x, y);
1895 :
1896 4411770 : case t_SER: {
1897 4411770 : vx = varn(x);
1898 4411770 : vy = varn(y);
1899 4411770 : if (vx != vy) {
1900 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1901 3675 : return mul_ser_scal(y, x);
1902 : }
1903 4408095 : lx = minss(lg(x), lg(y));
1904 4408095 : if (lx == 2) return zeroser(vx, valser(x)+valser(y));
1905 4408081 : av = avma; z = init_ser(lx, vx, valser(x)+valser(y));
1906 4408081 : x = ser2pol_i(x, lx);
1907 4408081 : y = ser2pol_i(y, lx);
1908 4408081 : y = RgXn_mul(x, y, lx-2);
1909 4408081 : return gerepilecopy(av, fill_ser(z,y));
1910 : }
1911 42 : case t_VEC:
1912 42 : if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
1913 : /* fall through, handle extended t_QFB */
1914 47804 : case t_QFB: return qfbcomp(x,y);
1915 6720472 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1916 4114579 : case t_MAT: return RgM_mul(x, y);
1917 :
1918 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1919 1631 : z = cgetg_copy(x, &l);
1920 1631 : if (l != lg(y)) break;
1921 17325 : for (i=1; i<l; i++)
1922 : {
1923 15694 : long yi = y[i];
1924 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1925 15694 : z[i] = x[yi];
1926 : }
1927 1631 : return z;
1928 :
1929 0 : default:
1930 0 : pari_err_TYPE2("*",x,y);
1931 : }
1932 : /* tx != ty */
1933 1759142603 : if (is_const_t(ty) && is_const_t(tx)) {
1934 1625605115 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1935 1625605115 : switch(tx) {
1936 1497270534 : case t_INT:
1937 : switch(ty)
1938 : {
1939 987265321 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1940 1322760 : case t_INTMOD:
1941 1322760 : z = cgetg(3, t_INTMOD);
1942 1322760 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1943 72286877 : case t_FRAC:
1944 72286877 : if (!signe(x)) return gen_0;
1945 55110010 : z=cgetg(3,t_FRAC);
1946 55110235 : p1 = gcdii(x,gel(y,2));
1947 55109297 : if (equali1(p1))
1948 : {
1949 28845681 : set_avma((pari_sp)z);
1950 28845672 : gel(z,2) = icopy(gel(y,2));
1951 28845758 : gel(z,1) = mulii(gel(y,1), x);
1952 : }
1953 : else
1954 : {
1955 26263954 : x = diviiexact(x,p1); tetpil = avma;
1956 26263305 : gel(z,2) = diviiexact(gel(y,2), p1);
1957 26263022 : gel(z,1) = mulii(gel(y,1), x);
1958 26263481 : fix_frac_if_int_GC(z,tetpil);
1959 : }
1960 55109809 : return z;
1961 433488111 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1962 4196299 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1963 1701 : case t_QUAD: return mulRq(x,y);
1964 1627719 : case t_FFELT: return FF_Z_mul(y,x);
1965 : }
1966 :
1967 : case t_REAL:
1968 120866840 : switch(ty)
1969 : {
1970 37415847 : case t_FRAC: return mulrfrac(x, y);
1971 83450951 : case t_COMPLEX: return mulRc(x, y);
1972 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1973 21 : default: pari_err_TYPE2("*",x,y);
1974 : }
1975 :
1976 8022 : case t_INTMOD:
1977 : switch(ty)
1978 : {
1979 7133 : case t_FRAC: { GEN X = gel(x,1);
1980 7133 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1981 7133 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1982 : }
1983 49 : case t_COMPLEX: return mulRc_direct(x,y);
1984 427 : case t_PADIC: { GEN X = gel(x,1);
1985 427 : z = cgetg(3, t_INTMOD);
1986 427 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1987 : }
1988 63 : case t_QUAD: return mulRq(x, y);
1989 350 : case t_FFELT:
1990 350 : if (!equalii(gel(x,1),FF_p_i(y)))
1991 0 : pari_err_OP("*",x,y);
1992 350 : return FF_Z_mul(y,gel(x,2));
1993 : }
1994 :
1995 : case t_FRAC:
1996 : switch(ty)
1997 : {
1998 5208322 : case t_COMPLEX: return mulRc(x, y);
1999 2581576 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
2000 343 : case t_QUAD: return mulRq(x, y);
2001 2051 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
2002 : }
2003 :
2004 : case t_FFELT:
2005 37 : pari_err_TYPE2("*",x,y);
2006 :
2007 21 : case t_COMPLEX:
2008 : switch(ty)
2009 : {
2010 14 : case t_PADIC:
2011 14 : return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
2012 7 : case t_QUAD:
2013 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
2014 7 : return mulqf(y, x, lx);
2015 : }
2016 :
2017 : case t_PADIC: /* ty == t_QUAD */
2018 28 : return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
2019 : }
2020 : }
2021 :
2022 135324739 : if (is_matvec_t(ty))
2023 : {
2024 8362556 : if (!is_matvec_t(tx))
2025 : {
2026 7413264 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
2027 7413271 : z = cgetg_copy(y, &ly);
2028 631997248 : for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
2029 7413036 : return z;
2030 : }
2031 949282 : switch(tx)
2032 : {
2033 122478 : case t_VEC:
2034 : switch(ty) {
2035 15540 : case t_COL: return RgV_RgC_mul(x,y);
2036 106938 : case t_MAT: return RgV_RgM_mul(x,y);
2037 : }
2038 0 : break;
2039 1687 : case t_COL:
2040 : switch(ty) {
2041 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2042 0 : case t_MAT: return RgC_RgM_mul(x,y);
2043 : }
2044 0 : break;
2045 825143 : case t_MAT:
2046 : switch(ty) {
2047 0 : case t_VEC: return RgM_RgV_mul(x,y);
2048 825143 : case t_COL: return RgM_RgC_mul(x,y);
2049 : }
2050 : }
2051 : }
2052 130386690 : if (is_matvec_t(tx))
2053 : {
2054 3058775 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2055 3059119 : z = cgetg_copy(x, &lx);
2056 10900791 : for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
2057 3058582 : return z;
2058 : }
2059 127328578 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2060 : /* tx < ty, !ismatvec(x and y) */
2061 :
2062 127328578 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2063 3404159 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2064 123924419 : if (is_scalar_t(tx)) {
2065 120633238 : if (tx == t_POLMOD) {
2066 3106889 : vx = varn(gel(x,1));
2067 3106889 : vy = gvar(y);
2068 3106889 : if (vx != vy) {
2069 2872550 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2070 23324 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2071 : }
2072 : /* error if ty == t_SER */
2073 234339 : av = avma; y = gmod(y, gel(x,1));
2074 234332 : return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2075 : }
2076 117526349 : return mul_scal(y, x, ty);
2077 : }
2078 :
2079 : /* x and y are not scalars, nor matvec */
2080 3291057 : vx = gvar(x);
2081 3291076 : vy = gvar(y);
2082 3291076 : if (vx != vy) /* x or y is treated as a scalar */
2083 2780236 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2084 2780236 : : mul_scal(y, x, ty);
2085 : /* vx = vy */
2086 1872181 : switch(tx)
2087 : {
2088 1872146 : case t_POL:
2089 : switch (ty)
2090 : {
2091 6741 : case t_SER:
2092 : {
2093 : long v;
2094 6741 : av = avma; v = RgX_valrem(x, &x);
2095 6741 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2096 6727 : v += valser(y); ly = lg(y);
2097 6727 : if (ly == 2) { set_avma(av); return zeroser(vx, v); }
2098 6496 : if (degpol(x))
2099 : {
2100 2030 : z = init_ser(ly, vx, v);
2101 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2102 2030 : return gerepilecopy(av, fill_ser(z, x));
2103 : }
2104 : /* take advantage of x = c*t^v */
2105 4466 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2106 4466 : setvalser(y, v); return y;
2107 : }
2108 :
2109 1865405 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2110 : }
2111 0 : break;
2112 :
2113 35 : case t_SER:
2114 : switch (ty)
2115 : {
2116 35 : case t_RFRAC:
2117 35 : av = avma;
2118 35 : return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2119 : }
2120 0 : break;
2121 : }
2122 0 : pari_err_TYPE2("*",x,y);
2123 : return NULL; /* LCOV_EXCL_LINE */
2124 : }
2125 :
2126 : /* return a nonnormalized result */
2127 : GEN
2128 112731 : sqr_ser_part(GEN x, long l1, long l2)
2129 : {
2130 : long i, j, l;
2131 : pari_sp av;
2132 : GEN Z, z, p1, p2;
2133 : long mi;
2134 112731 : if (l2 < l1) return zeroser(varn(x), 2*valser(x));
2135 112717 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2136 112717 : Z = cgetg(l2-l1+3, t_SER);
2137 112717 : Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
2138 112717 : z = Z + 2-l1;
2139 112717 : x += 2; mi = 0;
2140 425390 : for (i=0; i<l1; i++)
2141 : {
2142 312673 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2143 : }
2144 :
2145 720212 : for (i=l1; i<=l2; i++)
2146 : {
2147 607495 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2148 607495 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2149 1256134 : for (j=i-mi; j<=minss(l,mi); j++)
2150 648639 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2151 607495 : p1 = gshift(p1,1);
2152 607495 : if ((i&1) == 0 && p2[i>>1])
2153 93585 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2154 607495 : gel(z,i) = gerepileupto(av,p1);
2155 : }
2156 112717 : return Z;
2157 : }
2158 :
2159 : GEN
2160 1122028735 : gsqr(GEN x)
2161 : {
2162 : long i, lx;
2163 : pari_sp av, tetpil;
2164 : GEN z, p1, p2, p3, p4;
2165 :
2166 1122028735 : switch(typ(x))
2167 : {
2168 909745110 : case t_INT: return sqri(x);
2169 195041389 : case t_REAL: return sqrr(x);
2170 142489 : case t_INTMOD: { GEN X = gel(x,1);
2171 142489 : z = cgetg(3,t_INTMOD);
2172 142489 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
2173 142487 : gel(z,1) = icopy(X); return z;
2174 : }
2175 4140889 : case t_FRAC: return sqrfrac(x);
2176 :
2177 8085223 : case t_COMPLEX:
2178 8085223 : if (isintzero(gel(x,1))) {
2179 238806 : av = avma;
2180 238806 : return gerepileupto(av, gneg(gsqr(gel(x,2))));
2181 : }
2182 7846415 : z = cgetg(3,t_COMPLEX); av = avma;
2183 7846410 : p1 = gadd(gel(x,1),gel(x,2));
2184 7846225 : p2 = gsub(gel(x,1), gel(x,2));
2185 7846231 : p3 = gmul(gel(x,1),gel(x,2));
2186 7846334 : tetpil = avma;
2187 7846334 : gel(z,1) = gmul(p1,p2);
2188 7846379 : gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
2189 :
2190 4214 : case t_PADIC:
2191 4214 : z = cgetg(5,t_PADIC);
2192 4214 : i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
2193 4214 : if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
2194 4214 : z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
2195 4214 : gel(z,2) = icopy(gel(x,2));
2196 4214 : gel(z,3) = shifti(gel(x,3), i); av = avma;
2197 4214 : gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
2198 4214 : return z;
2199 :
2200 70 : case t_QUAD: z = cgetg(4,t_QUAD);
2201 70 : p1 = gel(x,1);
2202 70 : gel(z,1) = ZX_copy(p1); av = avma;
2203 70 : p2 = gsqr(gel(x,2));
2204 70 : p3 = gsqr(gel(x,3));
2205 70 : p4 = gmul(gneg_i(gel(p1,2)),p3);
2206 :
2207 70 : if (gequal0(gel(p1,3)))
2208 : {
2209 7 : tetpil = avma;
2210 7 : gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
2211 7 : av = avma;
2212 7 : p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
2213 7 : gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
2214 : }
2215 :
2216 63 : p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
2217 63 : tetpil = avma;
2218 63 : gel(z,2) = gadd(p2,p4);
2219 63 : gel(z,3) = gadd(p1,p3);
2220 63 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
2221 :
2222 484324 : case t_POLMOD:
2223 484324 : return sqr_polmod(gel(x,1), gel(x,2));
2224 :
2225 2957476 : case t_FFELT: return FF_sqr(x);
2226 :
2227 1327285 : case t_POL: return RgX_sqr(x);
2228 :
2229 35014 : case t_SER:
2230 35014 : lx = lg(x);
2231 35014 : if (ser_isexactzero(x)) {
2232 21 : GEN z = gcopy(x);
2233 21 : setvalser(z, 2*valser(x));
2234 21 : return z;
2235 : }
2236 34993 : if (lx < 40)
2237 34734 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2238 : else
2239 : {
2240 259 : pari_sp av = avma;
2241 259 : GEN z = init_ser(lx, varn(x), 2*valser(x));
2242 259 : x = ser2pol_i(x, lx);
2243 259 : x = RgXn_sqr(x, lx-2);
2244 259 : return gerepilecopy(av, fill_ser(z,x));
2245 : }
2246 :
2247 49 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2248 49 : gel(z,1) = gsqr(gel(x,1));
2249 49 : gel(z,2) = gsqr(gel(x,2)); return z;
2250 :
2251 1113 : case t_MAT: return RgM_sqr(x);
2252 28 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
2253 : /* fall through handle extended t_QFB */
2254 81919 : case t_QFB: return qfbsqr(x);
2255 658 : case t_VECSMALL:
2256 658 : z = cgetg_copy(x, &lx);
2257 16289 : for (i=1; i<lx; i++)
2258 : {
2259 15631 : long xi = x[i];
2260 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2261 15631 : z[i] = x[xi];
2262 : }
2263 658 : return z;
2264 : }
2265 7 : pari_err_TYPE2("*",x,x);
2266 : return NULL; /* LCOV_EXCL_LINE */
2267 : }
2268 :
2269 : /********************************************************************/
2270 : /** **/
2271 : /** DIVISION **/
2272 : /** **/
2273 : /********************************************************************/
2274 : static GEN
2275 13169 : div_rfrac_scal(GEN x, GEN y)
2276 : {
2277 13169 : pari_sp av = avma;
2278 13169 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2279 13169 : return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
2280 : }
2281 : static GEN
2282 37466 : div_scal_rfrac(GEN x, GEN y)
2283 : {
2284 37466 : GEN y1 = gel(y,1), y2 = gel(y,2);
2285 37466 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2286 : {
2287 14 : if (degpol(y1))
2288 : {
2289 14 : pari_sp av = avma;
2290 14 : GEN _1 = Rg_get_1(x);
2291 14 : if (transtype(_1)) y1 = gmul(y1, _1);
2292 14 : return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
2293 : }
2294 0 : y1 = gel(y1,2);
2295 : }
2296 37452 : return RgX_Rg_mul(y2, gdiv(x,y1));
2297 : }
2298 : static GEN
2299 1186663 : div_rfrac(GEN x, GEN y)
2300 1186663 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2301 :
2302 : /* x != 0 */
2303 : static GEN
2304 1338593 : div_ser_scal(GEN y, GEN x)
2305 : {
2306 : long i, l;
2307 : GEN z;
2308 1338593 : if (ser_isexactzero(y))
2309 : {
2310 28 : z = scalarser(lg(y) == 2? Rg_get_0(x): gdiv(gel(y,2), x), varn(y), 1);
2311 28 : setvalser(z, valser(y)); return z;
2312 : }
2313 1338565 : z = cgetg_copy(y, &l); z[1] = y[1];
2314 6202113 : for (i = 2; i < l; i++) gel(z,i) = gdiv(gel(y,i), x);
2315 1338565 : return normalizeser(z);
2316 : }
2317 : GEN
2318 7630 : ser_normalize(GEN x)
2319 : {
2320 7630 : long i, lx = lg(x);
2321 : GEN c, z;
2322 7630 : if (lx == 2) return x;
2323 7630 : c = gel(x,2); if (gequal1(c)) return x;
2324 7553 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2325 108528 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2326 7553 : return z;
2327 : }
2328 :
2329 : /* y != 0 */
2330 : static GEN
2331 6684328 : div_T_scal(GEN x, GEN y, long tx) {
2332 6684328 : switch(tx)
2333 : {
2334 5334646 : case t_POL: return RgX_Rg_div(x, y);
2335 1338586 : case t_SER: return div_ser_scal(x, y);
2336 11097 : case t_RFRAC: return div_rfrac_scal(x,y);
2337 : }
2338 0 : pari_err_TYPE2("/",x,y);
2339 : return NULL; /* LCOV_EXCL_LINE */
2340 : }
2341 :
2342 : static GEN
2343 9257463 : div_scal_pol(GEN x, GEN y) {
2344 9257463 : long ly = lg(y);
2345 : pari_sp av;
2346 : GEN _1;
2347 9257463 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2348 9179895 : if (isrationalzero(x)) return zeropol(varn(y));
2349 7130095 : av = avma;
2350 7130095 : _1 = Rg_get_1(x); if (transtype(_1)) y = gmul(y, _1);
2351 7130095 : return gerepileupto(av, gred_rfrac_simple(x,y));
2352 : }
2353 : static GEN
2354 18550 : div_scal_ser(GEN x, GEN y)
2355 : {
2356 18550 : pari_sp av = avma;
2357 18550 : GEN _1 = Rg_get_1(x);
2358 18550 : if (transtype(_1)) y = gmul(y, _1);
2359 18550 : return gerepileupto(av, gmul(x, ser_inv(y)));
2360 : }
2361 : static GEN
2362 9266775 : div_scal_T(GEN x, GEN y, long ty) {
2363 9266775 : switch(ty)
2364 : {
2365 9211522 : case t_POL: return div_scal_pol(x, y);
2366 18543 : case t_SER: return div_scal_ser(x, y);
2367 36710 : case t_RFRAC: return div_scal_rfrac(x, y);
2368 : }
2369 0 : pari_err_TYPE2("/",x,y);
2370 : return NULL; /* LCOV_EXCL_LINE */
2371 : }
2372 :
2373 : /* assume tx = ty = t_SER, same variable vx */
2374 : static GEN
2375 771496 : div_ser(GEN x, GEN y, long vx)
2376 : {
2377 771496 : long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
2378 771496 : GEN y0 = y, z;
2379 771496 : pari_sp av = avma;
2380 :
2381 771496 : if (!signe(y)) pari_err_INV("div_ser", y);
2382 771489 : if (ser_isexactzero(x))
2383 : {
2384 59892 : if (lx == 2) return zeroser(vx, v);
2385 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2386 147 : setvalser(z, v); return z;
2387 : }
2388 711597 : if (lx < ly) ly = lx;
2389 711597 : y = ser2pol_i_normalize(y, ly, &e);
2390 711597 : if (e)
2391 : {
2392 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2393 7 : ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
2394 : }
2395 711597 : z = init_ser(ly, vx, v);
2396 711597 : if (ly == 3)
2397 : {
2398 12333 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2399 12333 : if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
2400 12333 : return gerepileupto(av, z);
2401 : }
2402 699264 : x = ser2pol_i(x, ly);
2403 699264 : y = RgXn_div_i(x, y, ly-2);
2404 699264 : return gerepilecopy(av, fill_ser(z,y));
2405 : }
2406 : /* x,y compatible PADIC */
2407 : static GEN
2408 13210958 : divpp(GEN x, GEN y) {
2409 : pari_sp av;
2410 : long a, b;
2411 : GEN z, M;
2412 :
2413 13210958 : if (!signe(gel(y,4))) pari_err_INV("divpp",y);
2414 13210951 : if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
2415 13210622 : a = precp(x);
2416 13210622 : b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
2417 13210622 : z = cgetg(5, t_PADIC);
2418 13210238 : z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
2419 13210147 : gel(z,2) = icopy(gel(x,2));
2420 13209903 : gel(z,3) = icopy(M); av = avma;
2421 13209600 : gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
2422 13210450 : return z;
2423 : }
2424 : static GEN
2425 36736 : div_polmod_same(GEN T, GEN x, GEN y)
2426 : {
2427 36736 : long v = varn(T);
2428 36736 : GEN a, z = cgetg(3, t_POLMOD);
2429 36736 : gel(z,1) = RgX_copy(T);
2430 36736 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2431 11102 : a = gdiv(x, y);
2432 25634 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2433 119 : {
2434 119 : pari_sp av = avma;
2435 119 : a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
2436 : }
2437 25515 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2438 5271 : {
2439 5271 : pari_sp av = avma;
2440 5271 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2441 5271 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2442 5271 : a = gerepileupto(av, a);
2443 : }
2444 : else
2445 : {
2446 20244 : pari_sp av = avma;
2447 20244 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2448 20244 : a = gerepileupto(av, a);
2449 : }
2450 36736 : gel(z,2) = a; return z;
2451 : }
2452 : GEN
2453 432237464 : gdiv(GEN x, GEN y)
2454 : {
2455 432237464 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2456 : pari_sp av, tetpil;
2457 : GEN z, p1, p2;
2458 :
2459 432237464 : if (tx == ty) switch(tx)
2460 : {
2461 88643284 : case t_INT:
2462 88643284 : return Qdivii(x,y);
2463 :
2464 93476922 : case t_REAL: return divrr(x,y);
2465 25647 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2466 25647 : z = cgetg(3,t_INTMOD);
2467 25647 : if (X==Y || equalii(X,Y))
2468 25633 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2469 14 : gel(z,1) = gcdii(X,Y);
2470 14 : warn_coercion(X,Y,gel(z,1));
2471 14 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2472 14 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
2473 : }
2474 4040978 : case t_FRAC: {
2475 4040978 : GEN x1 = gel(x,1), x2 = gel(x,2);
2476 4040978 : GEN y1 = gel(y,1), y2 = gel(y,2);
2477 4040978 : z = cgetg(3, t_FRAC);
2478 4040979 : p1 = gcdii(x1, y1);
2479 4040977 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2480 4040977 : p1 = gcdii(x2, y2);
2481 4040977 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2482 4040972 : tetpil = avma;
2483 4040972 : gel(z,2) = mulii(x2,y1);
2484 4040972 : gel(z,1) = mulii(x1,y2);
2485 4040973 : normalize_frac(z);
2486 4040973 : fix_frac_if_int_GC(z,tetpil);
2487 4040980 : return z;
2488 : }
2489 9149130 : case t_COMPLEX:
2490 9149130 : if (isintzero(gel(y,1)))
2491 : {
2492 198002 : y = gel(y,2);
2493 198002 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2494 121093 : z = cgetg(3,t_COMPLEX);
2495 121093 : gel(z,1) = gdiv(gel(x,2), y);
2496 121093 : av = avma;
2497 121093 : gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
2498 121093 : return z;
2499 : }
2500 8951126 : av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
2501 8951095 : return gerepile(av, tetpil, gdiv(p2,p1));
2502 :
2503 587154 : case t_PADIC:
2504 587154 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
2505 587154 : return divpp(x, y);
2506 :
2507 238 : case t_QUAD:
2508 238 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2509 224 : av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
2510 224 : return gerepile(av, tetpil, gdiv(p2,p1));
2511 :
2512 243606 : case t_FFELT: return FF_div(x,y);
2513 :
2514 41111 : case t_POLMOD:
2515 41111 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2516 36736 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2517 : else {
2518 4375 : av = avma;
2519 4375 : z = gerepileupto(av, gmul(x, ginv(y)));
2520 : }
2521 41111 : return z;
2522 :
2523 18452705 : case t_POL:
2524 18452705 : vx = varn(x);
2525 18452705 : vy = varn(y);
2526 18452705 : if (vx != vy) {
2527 101962 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2528 45941 : else return div_scal_pol(x, y);
2529 : }
2530 18350743 : if (!signe(y)) pari_err_INV("gdiv",y);
2531 18350743 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2532 18228521 : av = avma;
2533 18228521 : return gerepileupto(av, gred_rfrac2(x,y));
2534 :
2535 771517 : case t_SER:
2536 771517 : vx = varn(x);
2537 771517 : vy = varn(y);
2538 771517 : if (vx != vy) {
2539 21 : if (varncmp(vx, vy) < 0)
2540 : {
2541 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2542 7 : return div_ser_scal(x, y);
2543 : }
2544 7 : return div_scal_ser(x, y);
2545 : }
2546 771496 : return div_ser(x, y, vx);
2547 1189351 : case t_RFRAC:
2548 1189351 : vx = varn(gel(x,2));
2549 1189351 : vy = varn(gel(y,2));
2550 1189351 : if (vx != vy) {
2551 2688 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2552 756 : else return div_scal_rfrac(x, y);
2553 : }
2554 1186663 : return div_rfrac(x,y);
2555 :
2556 21 : case t_VEC: /* handle extended t_QFB */
2557 21 : case t_QFB: av = avma; return gerepileupto(av, qfbcomp(x, ginv(y)));
2558 :
2559 28 : case t_MAT:
2560 28 : p1 = RgM_div(x,y);
2561 28 : if (!p1) pari_err_INV("gdiv",y);
2562 21 : return p1;
2563 :
2564 0 : default: pari_err_TYPE2("/",x,y);
2565 : }
2566 :
2567 215620824 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2568 : {
2569 3772884 : long s = signe(x);
2570 3772884 : if (!s) {
2571 516485 : if (gequal0(y)) pari_err_INV("gdiv",y);
2572 516485 : switch (ty)
2573 : {
2574 513321 : default: return gen_0;
2575 28 : case t_INTMOD:
2576 28 : z = cgetg(3,t_INTMOD);
2577 28 : gel(z,1) = icopy(gel(y,1));
2578 28 : gel(z,2) = gen_0; return z;
2579 3136 : case t_FFELT: return FF_zero(y);
2580 : }
2581 : }
2582 3256399 : if (is_pm1(x)) {
2583 1307543 : if (s > 0) return ginv(y);
2584 221003 : av = avma; return gerepileupto(av, ginv(gneg(y)));
2585 : }
2586 1948856 : switch(ty)
2587 : {
2588 638767 : case t_REAL: return divir(x,y);
2589 21 : case t_INTMOD:
2590 21 : z = cgetg(3, t_INTMOD);
2591 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2592 794523 : case t_FRAC:
2593 794523 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2594 794523 : if (equali1(p1))
2595 : {
2596 273744 : set_avma((pari_sp)z);
2597 273744 : gel(z,2) = icopy(gel(y,1));
2598 273744 : gel(z,1) = mulii(gel(y,2), x);
2599 273744 : normalize_frac(z);
2600 273744 : fix_frac_if_int(z);
2601 : }
2602 : else
2603 : {
2604 520779 : x = diviiexact(x,p1); tetpil = avma;
2605 520779 : gel(z,2) = diviiexact(gel(y,1), p1);
2606 520779 : gel(z,1) = mulii(gel(y,2), x);
2607 520779 : normalize_frac(z);
2608 520779 : fix_frac_if_int_GC(z,tetpil);
2609 : }
2610 794523 : return z;
2611 :
2612 469 : case t_FFELT: return Z_FF_div(x,y);
2613 513347 : case t_COMPLEX: return divRc(x,y);
2614 1505 : case t_PADIC: return divTp(x, y);
2615 231 : case t_QUAD:
2616 231 : av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
2617 231 : return gerepile(av, tetpil, gdiv(p2,p1));
2618 : }
2619 : }
2620 211847933 : if (gequal0(y))
2621 : {
2622 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2623 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2624 : }
2625 :
2626 211858176 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2627 : {
2628 172558913 : case t_REAL:
2629 172558913 : switch(ty)
2630 : {
2631 170198171 : case t_INT: return divri(x,y);
2632 522503 : case t_FRAC:
2633 522503 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2634 522503 : return gerepileuptoleaf(av, z);
2635 1838295 : case t_COMPLEX: return divRc(x, y);
2636 42 : case t_QUAD: return divfq(x, y, realprec(x));
2637 18 : default: pari_err_TYPE2("/",x,y);
2638 : }
2639 :
2640 280 : case t_INTMOD:
2641 : switch(ty)
2642 : {
2643 196 : case t_INT:
2644 196 : z = cgetg(3, t_INTMOD);
2645 196 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2646 28 : case t_FRAC: { GEN X = gel(x,1);
2647 28 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2648 28 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2649 : }
2650 7 : case t_FFELT:
2651 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2652 0 : pari_err_OP("/",x,y);
2653 7 : return Z_FF_div(gel(x,2),y);
2654 :
2655 0 : case t_COMPLEX:
2656 0 : av = avma;
2657 0 : return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
2658 :
2659 14 : case t_QUAD:
2660 14 : av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
2661 14 : return gerepile(av,tetpil, gdiv(p2,p1));
2662 :
2663 7 : case t_PADIC: { GEN X = gel(x,1);
2664 7 : z = cgetg(3, t_INTMOD);
2665 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2666 : }
2667 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2668 : }
2669 :
2670 : case t_FRAC:
2671 : switch(ty)
2672 : {
2673 2097442 : case t_INT: z = cgetg(3, t_FRAC);
2674 2097442 : p1 = gcdii(y,gel(x,1));
2675 2097442 : if (equali1(p1))
2676 : {
2677 827989 : set_avma((pari_sp)z); tetpil = 0;
2678 827989 : gel(z,1) = icopy(gel(x,1));
2679 : }
2680 : else
2681 : {
2682 1269453 : y = diviiexact(y,p1); tetpil = avma;
2683 1269453 : gel(z,1) = diviiexact(gel(x,1), p1);
2684 : }
2685 2097442 : gel(z,2) = mulii(gel(x,2),y);
2686 2097442 : normalize_frac(z);
2687 2097442 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2688 2097442 : return z;
2689 :
2690 59502 : case t_REAL:
2691 59502 : av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
2692 59502 : return gerepile(av, tetpil, divir(gel(x,1), p1));
2693 :
2694 7 : case t_INTMOD: { GEN Y = gel(y,1);
2695 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2696 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2697 : }
2698 :
2699 28 : case t_FFELT: av=avma;
2700 28 : return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2701 :
2702 861 : case t_COMPLEX: return divRc(x, y);
2703 :
2704 2141 : case t_PADIC:
2705 2141 : if (!signe(gel(x,1))) return gen_0;
2706 2141 : return divTp(x, y);
2707 :
2708 42 : case t_QUAD:
2709 42 : av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
2710 42 : return gerepile(av,tetpil,gdiv(p2,p1));
2711 : }
2712 :
2713 : case t_FFELT:
2714 161 : switch (ty)
2715 : {
2716 77 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2717 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2718 7 : case t_INTMOD:
2719 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2720 0 : pari_err_OP("/",x,y);
2721 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2722 49 : default:
2723 49 : pari_err_TYPE2("/",x,y);
2724 : }
2725 0 : break;
2726 :
2727 13452014 : case t_COMPLEX:
2728 : switch(ty)
2729 : {
2730 13452015 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2731 0 : case t_INTMOD: return mulRc_direct(ginv(y), x);
2732 0 : case t_PADIC:
2733 0 : return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
2734 0 : case t_QUAD:
2735 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2736 0 : return divfq(x, y, lx);
2737 : }
2738 :
2739 : case t_PADIC:
2740 : switch(ty)
2741 : {
2742 1204070 : case t_INT: case t_FRAC: { GEN p = gel(x,2);
2743 1203979 : return signe(gel(x,4))? divpT(x, y)
2744 2408049 : : zeropadic(p, valp(x) - Q_pval(y,p));
2745 : }
2746 7 : case t_INTMOD: { GEN Y = gel(y,1);
2747 7 : z = cgetg(3, t_INTMOD);
2748 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2749 : }
2750 14 : case t_COMPLEX: case t_QUAD:
2751 14 : av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
2752 14 : return gerepile(av,tetpil,gdiv(p1,p2));
2753 :
2754 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2755 : }
2756 :
2757 : case t_QUAD:
2758 : switch (ty)
2759 : {
2760 1190 : case t_INT: case t_INTMOD: case t_FRAC:
2761 1190 : z = cgetg(4,t_QUAD);
2762 1190 : gel(z,1) = ZX_copy(gel(x,1));
2763 1190 : gel(z,2) = gdiv(gel(x,2), y);
2764 1190 : gel(z,3) = gdiv(gel(x,3), y); return z;
2765 63 : case t_REAL: return divqf(x, y, realprec(y));
2766 14 : case t_PADIC: return divTp(x, y);
2767 0 : case t_COMPLEX:
2768 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2769 0 : return divqf(x, y, ly);
2770 : }
2771 : }
2772 22481178 : switch(ty) {
2773 581838 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2774 581838 : return gmul(x, ginv(y)); /* missing gerepile, for speed */
2775 42 : case t_MAT:
2776 42 : av = avma; p1 = RgM_inv(y);
2777 35 : if (!p1) pari_err_INV("gdiv",y);
2778 28 : return gerepileupto(av, gmul(x, p1));
2779 0 : case t_VEC: case t_COL:
2780 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2781 0 : pari_err_TYPE2("/",x,y);
2782 : }
2783 21899298 : switch(tx) {
2784 3838305 : case t_VEC: case t_COL: case t_MAT:
2785 3838305 : z = cgetg_copy(x, &lx);
2786 12961487 : for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
2787 3838251 : return z;
2788 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2789 0 : pari_err_TYPE2("/",x,y);
2790 : }
2791 :
2792 18060993 : vy = gvar(y);
2793 18061965 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2794 209668 : vx = varn(X);
2795 209668 : if (vx != vy) {
2796 209094 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2797 208765 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2798 : }
2799 : /* y is POL, SER or RFRAC */
2800 574 : av = avma;
2801 574 : switch(ty)
2802 : {
2803 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2804 574 : default: y = ginvmod(gmod(y,X), X);
2805 : }
2806 567 : return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
2807 : }
2808 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2809 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2810 : * variable NO_VARIABLE, then the operation is incorrect */
2811 17852297 : vx = gvar(x);
2812 17852290 : if (vx != vy) { /* includes cases where one is scalar */
2813 15950778 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2814 9266449 : else return div_scal_T(x, y, ty);
2815 : }
2816 1901512 : switch(tx)
2817 : {
2818 1046164 : case t_POL:
2819 : switch(ty)
2820 : {
2821 28 : case t_SER:
2822 : {
2823 28 : GEN y0 = y;
2824 : long v;
2825 28 : av = avma; v = RgX_valrem(x, &x);
2826 28 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2827 14 : v -= valser(y); ly = lg(y); /* > 2 */
2828 14 : y = ser2pol_i_normalize(y, ly, &i);
2829 14 : if (i)
2830 : {
2831 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2832 7 : ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
2833 : }
2834 7 : z = init_ser(ly, vx, v);
2835 7 : return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, ly-2)));
2836 : }
2837 :
2838 1046136 : case t_RFRAC:
2839 : {
2840 1046136 : GEN y1 = gel(y,1), y2 = gel(y,2);
2841 1046136 : if (typ(y1) == t_POL && varn(y1) == vx)
2842 1171 : return mul_rfrac_scal(y2, y1, x);
2843 1044965 : av = avma;
2844 1044965 : return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2845 : }
2846 : }
2847 0 : break;
2848 :
2849 584324 : case t_SER:
2850 : switch(ty)
2851 : {
2852 584310 : case t_POL:
2853 : {
2854 584310 : long v = valser(x);
2855 584310 : lx = lg(x);
2856 584310 : if (lx == 2) return zeroser(vx, v - RgX_val(y));
2857 584205 : av = avma;
2858 584205 : x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
2859 584205 : z = init_ser(lx, vx, v);
2860 584205 : if (!signe(x)) setsigne(z,0);
2861 584205 : return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
2862 : }
2863 14 : case t_RFRAC:
2864 14 : av = avma;
2865 14 : return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2866 : }
2867 0 : break;
2868 :
2869 271003 : case t_RFRAC:
2870 : switch(ty)
2871 : {
2872 271003 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2873 0 : case t_SER:
2874 0 : av = avma;
2875 0 : return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2876 : }
2877 0 : break;
2878 : }
2879 21 : pari_err_TYPE2("/",x,y);
2880 : return NULL; /* LCOV_EXCL_LINE */
2881 : }
2882 :
2883 : /********************************************************************/
2884 : /** **/
2885 : /** SIMPLE MULTIPLICATION **/
2886 : /** **/
2887 : /********************************************************************/
2888 : GEN
2889 229467168 : gmulsg(long s, GEN y)
2890 : {
2891 : long ly, i;
2892 : pari_sp av;
2893 : GEN z;
2894 :
2895 229467168 : switch(typ(y))
2896 : {
2897 136803119 : case t_INT: return mulsi(s,y);
2898 69360886 : case t_REAL: return s? mulsr(s,y): gen_0; /* gmul semantic */
2899 174118 : case t_INTMOD: { GEN p = gel(y,1);
2900 174118 : z = cgetg(3,t_INTMOD);
2901 174119 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
2902 174121 : gel(z,1) = icopy(p); return z;
2903 : }
2904 548202 : case t_FFELT: return FF_Z_mul(y,stoi(s));
2905 6204912 : case t_FRAC:
2906 6204912 : if (!s) return gen_0;
2907 6190268 : z = cgetg(3,t_FRAC);
2908 6190267 : i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
2909 6190267 : if (i == 1)
2910 : {
2911 2267707 : gel(z,2) = icopy(gel(y,2));
2912 2267707 : gel(z,1) = mulis(gel(y,1), s);
2913 : }
2914 : else
2915 : {
2916 3922560 : gel(z,2) = diviuexact(gel(y,2), (ulong)i);
2917 3922560 : gel(z,1) = mulis(gel(y,1), s/i);
2918 3922560 : fix_frac_if_int(z);
2919 : }
2920 6190267 : return z;
2921 :
2922 14198622 : case t_COMPLEX:
2923 14198622 : if (!s) return gen_0;
2924 13072818 : z = cgetg(3, t_COMPLEX);
2925 13072794 : gel(z,1) = gmulsg(s,gel(y,1));
2926 13071306 : gel(z,2) = gmulsg(s,gel(y,2)); return z;
2927 :
2928 1435 : case t_PADIC:
2929 1435 : if (!s) return gen_0;
2930 1435 : av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
2931 :
2932 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2933 7 : gel(z,1) = ZX_copy(gel(y,1));
2934 7 : gel(z,2) = gmulsg(s,gel(y,2));
2935 7 : gel(z,3) = gmulsg(s,gel(y,3)); return z;
2936 :
2937 205359 : case t_POLMOD:
2938 205359 : retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
2939 :
2940 746644 : case t_POL:
2941 746644 : if (!signe(y)) return RgX_copy(y);
2942 730327 : if (!s) return scalarpol(Rg_get_0(y), varn(y));
2943 727114 : z = cgetg_copy(y, &ly); z[1]=y[1];
2944 3084705 : for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2945 727113 : return normalizepol_lg(z, ly);
2946 :
2947 182 : case t_SER:
2948 182 : if (ser_isexactzero(y)) return gcopy(y);
2949 182 : if (!s) return Rg_get_0(y);
2950 182 : z = cgetg_copy(y, &ly); z[1]=y[1];
2951 3864 : for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2952 182 : return normalizeser(z);
2953 :
2954 0 : case t_RFRAC:
2955 0 : if (!s) return zeropol(varn(gel(y,2)));
2956 0 : if (s == 1) return gcopy(y);
2957 0 : if (s == -1) return gneg(y);
2958 0 : return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
2959 :
2960 1233392 : case t_VEC: case t_COL: case t_MAT:
2961 1233392 : z = cgetg_copy(y, &ly);
2962 3882086 : for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2963 1233390 : return z;
2964 : }
2965 0 : pari_err_TYPE("gmulsg",y);
2966 : return NULL; /* LCOV_EXCL_LINE */
2967 : }
2968 :
2969 : GEN
2970 124319643 : gmulug(ulong s, GEN y)
2971 : {
2972 : long ly, i;
2973 : pari_sp av;
2974 : GEN z;
2975 :
2976 124319643 : switch(typ(y))
2977 : {
2978 122204201 : case t_INT: return mului(s,y);
2979 1059371 : case t_REAL: return s? mulur(s,y): gen_0; /* gmul semantic */
2980 364 : case t_INTMOD: { GEN p = gel(y,1);
2981 364 : z = cgetg(3,t_INTMOD);
2982 364 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(mului(s,gel(y,2)), p));
2983 364 : gel(z,1) = icopy(p); return z;
2984 : }
2985 413 : case t_FFELT: return FF_Z_mul(y,utoi(s));
2986 979912 : case t_FRAC:
2987 979912 : if (!s) return gen_0;
2988 979898 : z = cgetg(3,t_FRAC);
2989 979898 : i = ugcd(s, umodiu(gel(y,2), s));
2990 979898 : if (i == 1)
2991 : {
2992 821635 : gel(z,2) = icopy(gel(y,2));
2993 821635 : gel(z,1) = muliu(gel(y,1), s);
2994 : }
2995 : else
2996 : {
2997 158263 : gel(z,2) = diviuexact(gel(y,2), i);
2998 158263 : gel(z,1) = muliu(gel(y,1), s/i);
2999 158263 : fix_frac_if_int(z);
3000 : }
3001 979898 : return z;
3002 :
3003 30765 : case t_COMPLEX:
3004 30765 : if (!s) return gen_0;
3005 30765 : z = cgetg(3, t_COMPLEX);
3006 30765 : gel(z,1) = gmulug(s,gel(y,1));
3007 30765 : gel(z,2) = gmulug(s,gel(y,2)); return z;
3008 :
3009 7342 : case t_PADIC:
3010 7342 : if (!s) return gen_0;
3011 7342 : av = avma; return gerepileupto(av, mulpp(cvtop2(utoi(s),y), y));
3012 :
3013 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3014 0 : gel(z,1) = ZX_copy(gel(y,1));
3015 0 : gel(z,2) = gmulug(s,gel(y,2));
3016 0 : gel(z,3) = gmulug(s,gel(y,3)); return z;
3017 :
3018 6783 : case t_POLMOD:
3019 6783 : retmkpolmod(gmulug(s,gel(y,2)), RgX_copy(gel(y,1)));
3020 :
3021 16226 : case t_POL:
3022 16226 : if (!signe(y)) return RgX_copy(y);
3023 15449 : if (!s) return scalarpol(Rg_get_0(y), varn(y));
3024 15449 : z = cgetg_copy(y, &ly); z[1]=y[1];
3025 52843 : for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
3026 15449 : return normalizepol_lg(z, ly);
3027 :
3028 0 : case t_SER:
3029 0 : if (ser_isexactzero(y)) return gcopy(y);
3030 0 : if (!s) return Rg_get_0(y);
3031 0 : z = cgetg_copy(y, &ly); z[1]=y[1];
3032 0 : for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
3033 0 : return normalizeser(z);
3034 :
3035 0 : case t_RFRAC:
3036 0 : if (!s) return zeropol(varn(gel(y,2)));
3037 0 : if (s == 1) return gcopy(y);
3038 0 : return mul_rfrac_scal(gel(y,1), gel(y,2), utoi(s));
3039 :
3040 14266 : case t_VEC: case t_COL: case t_MAT:
3041 14266 : z = cgetg_copy(y, &ly);
3042 74284 : for (i=1; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
3043 14266 : return z;
3044 : }
3045 0 : pari_err_TYPE("gmulsg",y);
3046 : return NULL; /* LCOV_EXCL_LINE */
3047 : }
3048 :
3049 : /********************************************************************/
3050 : /** **/
3051 : /** SIMPLE DIVISION **/
3052 : /** **/
3053 : /********************************************************************/
3054 :
3055 : GEN
3056 12720394 : gdivgs(GEN x, long s)
3057 : {
3058 12720394 : long tx = typ(x), lx, i;
3059 : pari_sp av;
3060 : GEN z;
3061 :
3062 12720394 : if (!s)
3063 : {
3064 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3065 0 : pari_err_INV("gdivgs",gen_0);
3066 : }
3067 12720398 : switch(tx)
3068 : {
3069 1561293 : case t_INT: return Qdivis(x, s);
3070 8377105 : case t_REAL: return divrs(x,s);
3071 :
3072 357 : case t_INTMOD:
3073 357 : z = cgetg(3, t_INTMOD);
3074 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3075 :
3076 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3077 :
3078 552004 : case t_FRAC: z = cgetg(3, t_FRAC);
3079 552004 : i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
3080 552004 : if (i == 1)
3081 : {
3082 409183 : gel(z,2) = mulsi(s, gel(x,2));
3083 409183 : gel(z,1) = icopy(gel(x,1));
3084 : }
3085 : else
3086 : {
3087 142821 : gel(z,2) = mulsi(s/i, gel(x,2));
3088 142821 : gel(z,1) = divis(gel(x,1), i);
3089 : }
3090 552004 : normalize_frac(z);
3091 552004 : fix_frac_if_int(z); return z;
3092 :
3093 1781962 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3094 1781962 : gel(z,1) = gdivgs(gel(x,1),s);
3095 1781960 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3096 :
3097 133 : case t_PADIC: /* divpT */
3098 : {
3099 133 : GEN p = gel(x,2);
3100 133 : if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
3101 133 : av = avma;
3102 133 : return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
3103 : }
3104 :
3105 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3106 28 : gel(z,1) = ZX_copy(gel(x,1));
3107 28 : gel(z,2) = gdivgs(gel(x,2),s);
3108 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3109 :
3110 37828 : case t_POLMOD:
3111 37828 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3112 :
3113 91 : case t_RFRAC:
3114 91 : if (s == 1) return gcopy(x);
3115 84 : else if (s == -1) return gneg(x);
3116 84 : return div_rfrac_scal(x, stoi(s));
3117 :
3118 37625 : case t_POL: case t_SER:
3119 37625 : z = cgetg_copy(x, &lx); z[1] = x[1];
3120 157454 : for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
3121 37625 : return z;
3122 371237 : case t_VEC: case t_COL: case t_MAT:
3123 371237 : z = cgetg_copy(x, &lx);
3124 822265 : for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
3125 371237 : return z;
3126 :
3127 : }
3128 0 : pari_err_TYPE2("/",x, stoi(s));
3129 : return NULL; /* LCOV_EXCL_LINE */
3130 : }
3131 :
3132 : GEN
3133 54707126 : gdivgu(GEN x, ulong s)
3134 : {
3135 54707126 : long tx = typ(x), lx, i;
3136 : pari_sp av;
3137 : GEN z;
3138 :
3139 54707126 : if (!s)
3140 : {
3141 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3142 0 : pari_err_INV("gdivgu",gen_0);
3143 : }
3144 54707071 : switch(tx)
3145 : {
3146 17910832 : case t_INT: return Qdiviu(x, s);
3147 12954230 : case t_REAL: return divru(x,s);
3148 :
3149 210315 : case t_INTMOD:
3150 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3151 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3152 :
3153 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3154 :
3155 641842 : case t_FRAC: z = cgetg(3, t_FRAC);
3156 641842 : i = ugcd(s, umodiu(gel(x,1), s));
3157 641842 : if (i == 1)
3158 : {
3159 461994 : gel(z,2) = mului(s, gel(x,2));
3160 461994 : gel(z,1) = icopy(gel(x,1));
3161 : }
3162 : else
3163 : {
3164 179848 : gel(z,2) = mului(s/i, gel(x,2));
3165 179848 : gel(z,1) = divis(gel(x,1), i);
3166 : }
3167 641842 : normalize_frac(z);
3168 641842 : fix_frac_if_int(z); return z;
3169 :
3170 11383791 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3171 11383791 : gel(z,1) = gdivgu(gel(x,1),s);
3172 11383791 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3173 :
3174 11551118 : case t_PADIC: /* divpT */
3175 : {
3176 11551118 : GEN p = gel(x,2);
3177 11551118 : if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
3178 11415592 : av = avma;
3179 11415592 : return gerepileupto(av, divpp(x, cvtop2(utoi(s),x)));
3180 : }
3181 :
3182 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3183 0 : gel(z,1) = ZX_copy(gel(x,1));
3184 0 : gel(z,2) = gdivgu(gel(x,2),s);
3185 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3186 :
3187 1456 : case t_POLMOD:
3188 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3189 :
3190 56 : case t_RFRAC:
3191 56 : if (s == 1) return gcopy(x);
3192 56 : return div_rfrac_scal(x, utoi(s));
3193 :
3194 52794 : case t_POL: case t_SER:
3195 52794 : z = cgetg_copy(x, &lx); z[1] = x[1];
3196 222138 : for (i=2; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
3197 52794 : return z;
3198 329 : case t_VEC: case t_COL: case t_MAT:
3199 329 : z = cgetg_copy(x, &lx);
3200 1148 : for (i=1; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
3201 329 : return z;
3202 : }
3203 0 : pari_err_TYPE2("/",x, utoi(s));
3204 : return NULL; /* LCOV_EXCL_LINE */
3205 : }
3206 :
3207 : /* x / (i*(i+1)) */
3208 : GEN
3209 223274041 : divrunextu(GEN x, ulong i)
3210 : {
3211 223274041 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3212 0 : return divri(x, muluu(i , i+1));
3213 : else
3214 223274041 : return divru(x, i*(i+1));
3215 : }
3216 : /* x / (i*(i+1)) */
3217 : GEN
3218 813628 : gdivgunextu(GEN x, ulong i)
3219 : {
3220 813628 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3221 0 : return gdivgu(x, i*(i+1));
3222 : else
3223 813628 : return gdiv(x, muluu(i, i+1));
3224 : }
3225 :
3226 : /* True shift (exact multiplication by 2^n) */
3227 : GEN
3228 124630614 : gmul2n(GEN x, long n)
3229 : {
3230 : long lx, i, k, l;
3231 : GEN z, a, b;
3232 :
3233 124630614 : switch(typ(x))
3234 : {
3235 31203238 : case t_INT:
3236 31203238 : if (n>=0) return shifti(x,n);
3237 5033952 : if (!signe(x)) return gen_0;
3238 3340590 : l = vali(x); n = -n;
3239 3340653 : if (n<=l) return shifti(x,-n);
3240 397681 : z = cgetg(3,t_FRAC);
3241 397681 : gel(z,1) = shifti(x,-l);
3242 397681 : gel(z,2) = int2n(n-l); return z;
3243 :
3244 62756806 : case t_REAL:
3245 62756806 : return shiftr(x,n);
3246 :
3247 180971 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3248 180971 : z = cgetg(3,t_INTMOD);
3249 180971 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3250 82810 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
3251 82806 : gel(z,1) = icopy(b); return z;
3252 :
3253 217412 : case t_FFELT: return FF_mul2n(x,n);
3254 :
3255 1268222 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3256 1268222 : l = vali(a);
3257 1268222 : k = vali(b);
3258 1268222 : if (n+l >= k)
3259 : {
3260 395165 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3261 302676 : l = n-k; k = -k;
3262 : }
3263 : else
3264 : {
3265 873057 : k = -(l+n); l = -l;
3266 : }
3267 1175733 : z = cgetg(3,t_FRAC);
3268 1175733 : gel(z,1) = shifti(a,l);
3269 1175733 : gel(z,2) = shifti(b,k); return z;
3270 :
3271 10721603 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3272 10721605 : gel(z,1) = gmul2n(gel(x,1),n);
3273 10721601 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3274 :
3275 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3276 105 : gel(z,1) = ZX_copy(gel(x,1));
3277 105 : gel(z,2) = gmul2n(gel(x,2),n);
3278 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3279 :
3280 193921 : case t_POLMOD:
3281 193921 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3282 :
3283 1717796 : case t_POL:
3284 1717796 : z = cgetg_copy(x, &lx); z[1] = x[1];
3285 9105548 : for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3286 1717825 : return normalizepol_lg(z, lx); /* needed if char = 2 */
3287 102833 : case t_SER:
3288 102833 : if (ser_isexactzero(x)) return gcopy(x);
3289 102805 : z = cgetg_copy(x, &lx); z[1] = x[1];
3290 634228 : for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3291 102805 : return normalizeser(z); /* needed if char = 2 */
3292 16264399 : case t_VEC: case t_COL: case t_MAT:
3293 16264399 : z = cgetg_copy(x, &lx);
3294 59925980 : for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3295 16264389 : return z;
3296 :
3297 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3298 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3299 :
3300 5320 : case t_PADIC: /* int2n wrong if n < 0 */
3301 5320 : return gmul(gmul2n(gen_1,n),x);
3302 : }
3303 0 : pari_err_TYPE("gmul2n",x);
3304 : return NULL; /* LCOV_EXCL_LINE */
3305 : }
3306 :
3307 : /*******************************************************************/
3308 : /* */
3309 : /* INVERSE */
3310 : /* */
3311 : /*******************************************************************/
3312 : static GEN
3313 209275 : inv_polmod(GEN T, GEN x)
3314 : {
3315 209275 : GEN z = cgetg(3,t_POLMOD), a;
3316 209275 : gel(z,1) = RgX_copy(T);
3317 209273 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3318 84640 : a = ginv(x);
3319 : else
3320 : {
3321 124633 : if (lg(T) == 5) /* quadratic fields */
3322 13055 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3323 : else
3324 111578 : a = RgXQ_inv(x, T);
3325 : }
3326 209274 : gel(z,2) = a; return z;
3327 : }
3328 :
3329 : GEN
3330 36427734 : ginv(GEN x)
3331 : {
3332 : long s;
3333 : pari_sp av, tetpil;
3334 : GEN z, y, p1, p2;
3335 :
3336 36427734 : switch(typ(x))
3337 : {
3338 9498327 : case t_INT:
3339 9498327 : if (is_pm1(x)) return icopy(x);
3340 7894944 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3341 7894930 : z = cgetg(3,t_FRAC);
3342 7900353 : gel(z,1) = s<0? gen_m1: gen_1;
3343 7900353 : gel(z,2) = absi(x); return z;
3344 :
3345 10710766 : case t_REAL: return invr(x);
3346 :
3347 12964 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3348 12964 : gel(z,1) = icopy(gel(x,1));
3349 12964 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3350 :
3351 551982 : case t_FRAC: {
3352 551982 : GEN a = gel(x,1), b = gel(x,2);
3353 551982 : s = signe(a);
3354 551982 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3355 248264 : z = cgetg(3,t_FRAC);
3356 248264 : gel(z,1) = icopy(b);
3357 248264 : gel(z,2) = icopy(a);
3358 248264 : normalize_frac(z); return z;
3359 : }
3360 9582987 : case t_COMPLEX:
3361 9582987 : av=avma;
3362 9582987 : p1=cxnorm(x);
3363 9582217 : p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
3364 9582860 : tetpil=avma;
3365 9582860 : return gerepile(av,tetpil,divcR(p2,p1));
3366 :
3367 273 : case t_QUAD:
3368 273 : av=avma; p1=quadnorm(x); p2=conj_i(x); tetpil=avma;
3369 273 : return gerepile(av,tetpil,gdiv(p2,p1));
3370 :
3371 358085 : case t_PADIC: z = cgetg(5,t_PADIC);
3372 358085 : if (!signe(gel(x,4))) pari_err_INV("ginv",x);
3373 358078 : z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
3374 358078 : gel(z,2) = icopy(gel(x,2));
3375 358078 : gel(z,3) = icopy(gel(x,3));
3376 358078 : gel(z,4) = Zp_inv(gel(x,4),gel(z,2),precp(x)); return z;
3377 :
3378 209275 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3379 14215 : case t_FFELT: return FF_inv(x);
3380 5322709 : case t_POL: return gred_rfrac_simple(gen_1,x);
3381 26327 : case t_SER: return ser_inv(x);
3382 2926 : case t_RFRAC:
3383 : {
3384 2926 : GEN n = gel(x,1), d = gel(x,2);
3385 2926 : pari_sp av = avma, ltop;
3386 2926 : if (gequal0(n)) pari_err_INV("ginv",x);
3387 :
3388 2926 : n = simplify_shallow(n);
3389 2926 : if (typ(n) != t_POL || varn(n) != varn(d))
3390 : {
3391 2926 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3392 679 : ltop = avma;
3393 679 : z = RgX_Rg_div(d,n);
3394 : } else {
3395 0 : ltop = avma;
3396 0 : z = cgetg(3,t_RFRAC);
3397 0 : gel(z,1) = RgX_copy(d);
3398 0 : gel(z,2) = RgX_copy(n);
3399 : }
3400 679 : stackdummy(av, ltop);
3401 679 : return z;
3402 : }
3403 :
3404 21 : case t_VEC: if (!is_ext_qfr(x)) break;
3405 : case t_QFB:
3406 99521 : return qfbpow(x, gen_m1);
3407 38820 : case t_MAT:
3408 38820 : y = RgM_inv(x);
3409 38813 : if (!y) pari_err_INV("ginv",x);
3410 38743 : return y;
3411 28 : case t_VECSMALL:
3412 : {
3413 28 : long i, lx = lg(x)-1;
3414 28 : y = zero_zv(lx);
3415 112 : for (i=1; i<=lx; i++)
3416 : {
3417 84 : long xi = x[i];
3418 84 : if (xi<1 || xi>lx || y[xi])
3419 0 : pari_err_TYPE("ginv [not a permutation]", x);
3420 84 : y[xi] = i;
3421 : }
3422 28 : return y;
3423 : }
3424 : }
3425 6 : pari_err_TYPE("inverse",x);
3426 : return NULL; /* LCOV_EXCL_LINE */
3427 : }
|
