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 261167112 : addRc(GEN x, GEN y) {
61 261167112 : GEN z = cgetg(3,t_COMPLEX);
62 261169044 : gel(z,1) = gadd(x,gel(y,1));
63 261110348 : gel(z,2) = gcopy(gel(y,2)); return z;
64 : }
65 : static GEN
66 310080619 : mulRc(GEN x, GEN y) {
67 310080619 : GEN z = cgetg(3,t_COMPLEX);
68 310001322 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
69 309826225 : 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 2518162 : divRc(GEN x, GEN y) {
80 2518162 : GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
81 2518155 : GEN z = cgetg(3,t_COMPLEX);
82 2518158 : gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
83 2518145 : gel(z,2) = gmul(mt, gel(y,2));
84 2518141 : return z;
85 : }
86 : static GEN
87 19359055 : divcR(GEN x, GEN y) {
88 19359055 : GEN z = cgetg(3,t_COMPLEX);
89 19358997 : gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
90 19357995 : 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 30207420 : mulrfrac(GEN x, GEN y)
114 : {
115 : pari_sp av;
116 30207420 : GEN z, a = gel(y,1), b = gel(y,2);
117 30207420 : if (is_pm1(a)) /* frequent special case */
118 : {
119 5468388 : z = divri(x, b); if (signe(a) < 0) togglesign(z);
120 5468279 : return z;
121 : }
122 24739015 : av = avma;
123 24739015 : 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 6502491 : mulTp(GEN x, GEN y) { pari_sp av = avma;
151 6502491 : return gerepileupto(av, mulpp(cvtop2(x,y), y));
152 : }
153 : /* y PADIC, non zero x / y by converting x to padic */
154 : static GEN
155 3647 : divTp(GEN x, GEN y) { pari_sp av = avma;
156 3647 : return gerepileupto(av, divpp(cvtop2(x,y), y));
157 : }
158 : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
159 : * converted to O(p^e) and division by 0 */
160 : static GEN
161 1327543 : divpT(GEN x, GEN y) { pari_sp av = avma;
162 1327543 : 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 3634917 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
169 3634917 : if (lgefint(X) == 3) {
170 3367319 : ulong u = Fl_add(itou(x),itou(y), X[2]);
171 3367319 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
172 : }
173 : else {
174 267598 : GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
175 267598 : gel(z,2) = gerepileuptoint((pari_sp)z, u);
176 : }
177 3634917 : gel(z,1) = icopy(X); return z;
178 : }
179 : static GEN
180 1158469 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
181 1158469 : if (lgefint(X) == 3) {
182 784435 : ulong u = Fl_sub(itou(x),itou(y), X[2]);
183 784435 : 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 1158469 : gel(z,1) = icopy(X); return z;
190 : }
191 : /* cf add_intmod_same */
192 : static GEN
193 2936858 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
194 2936858 : if (lgefint(X) == 3) {
195 2795225 : ulong u = Fl_mul(itou(x),itou(y), X[2]);
196 2795225 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
197 : }
198 : else
199 141633 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
200 2936869 : gel(z,1) = icopy(X); return z;
201 : }
202 : /* cf add_intmod_same */
203 : static GEN
204 341858 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
205 : {
206 341858 : if (lgefint(X) == 3) {
207 319265 : ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
208 319258 : set_avma((pari_sp)z); gel(z,2) = utoi(u);
209 : }
210 : else
211 22593 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
212 341852 : 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 9891058 : rfrac_denom_mul_scal(GEN d, GEN y)
225 : {
226 9891058 : GEN D = RgX_Rg_mul(d, y);
227 9891058 : 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 9891058 : return D;
233 : }
234 :
235 : static int
236 57848015 : Leading_is_neg(GEN x)
237 : {
238 122327546 : while (typ(x) == t_POL) x = leading_coeff(x);
239 57848015 : return is_real_t(typ(x))? gsigne(x) < 0: 0;
240 : }
241 :
242 : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
243 : GEN
244 57852140 : gred_rfrac_simple(GEN n, GEN d)
245 : {
246 : GEN c, cn, cd, z;
247 57852140 : long dd = degpol(d);
248 :
249 57852140 : if (dd <= 0)
250 : {
251 4125 : if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
252 4125 : n = gdiv(n, gel(d,2));
253 4125 : if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
254 4125 : return n;
255 : }
256 57848015 : if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
257 57848015 : cd = content(d);
258 59719341 : while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
259 57848015 : cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
260 57848015 : if (!gequal1(cd)) {
261 6582193 : d = RgX_Rg_div(d,cd);
262 6582193 : if (!gequal1(cn))
263 : {
264 1309822 : if (gequal0(cn)) {
265 49 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
266 0 : n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
267 0 : c = gen_1;
268 : } else {
269 1309773 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
270 1309773 : c = gdiv(cn,cd);
271 : }
272 : }
273 : else
274 5272371 : c = ginv(cd);
275 : } else {
276 51265822 : if (!gequal1(cn))
277 : {
278 3279545 : if (gequal0(cn)) {
279 1484 : if (isexactzero(cn)) return scalarpol(cn, varn(d));
280 21 : c = gen_1;
281 : } else {
282 3278061 : n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
283 3278061 : c = cn;
284 : }
285 : } else {
286 47986277 : GEN y = cgetg(3,t_RFRAC);
287 47986277 : gel(y,1) = gcopy(n);
288 47986277 : gel(y,2) = RgX_copy(d); return y;
289 : }
290 : }
291 :
292 9860226 : if (typ(c) == t_POL)
293 : {
294 904069 : z = c;
295 943213 : do { z = content(z); } while (typ(z) == t_POL);
296 904069 : cd = denom_i(z);
297 904069 : cn = gmul(c, cd);
298 : }
299 : else
300 : {
301 8956157 : cn = numer_i(c);
302 8956157 : cd = denom_i(c);
303 : }
304 9860226 : z = cgetg(3,t_RFRAC);
305 9860226 : gel(z,1) = gmul(n, cn);
306 9860226 : gel(z,2) = d = rfrac_denom_mul_scal(d, cd);
307 : /* can happen: Pol(O(17^-1)) / Pol([Mod(9,23), O(23^-3)]) */
308 9860226 : if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
309 9860219 : return z;
310 : }
311 :
312 : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
313 : static GEN
314 248283 : fix_rfrac(GEN x, long d)
315 : {
316 : GEN z, N, D;
317 248283 : if (!d || typ(x) == t_POL) return x;
318 171682 : z = cgetg(3, t_RFRAC);
319 171682 : N = gel(x,1);
320 171682 : D = gel(x,2);
321 171682 : if (d > 0) {
322 8694 : gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
323 165893 : : monomialcopy(N,d,varn(D));
324 165844 : gel(z, 2) = RgX_copy(D);
325 : } else {
326 5838 : gel(z, 1) = gcopy(N);
327 5838 : gel(z, 2) = RgX_shift(D, -d);
328 : }
329 171682 : return z;
330 : }
331 :
332 : /* assume d != 0 */
333 : static GEN
334 44775710 : gred_rfrac2(GEN n, GEN d)
335 : {
336 : GEN y, z;
337 : long v, vd, vn;
338 :
339 44775710 : n = simplify_shallow(n);
340 44775710 : if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
341 37944074 : d = simplify_shallow(d);
342 37944074 : if (typ(d) != t_POL) return gdiv(n,d);
343 36764257 : vd = varn(d);
344 36764257 : if (typ(n) != t_POL)
345 : {
346 20625719 : if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
347 20624312 : if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
348 0 : pari_err_BUG("gred_rfrac2 [incompatible variables]");
349 : }
350 16138538 : vn = varn(n);
351 16138538 : if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
352 16000892 : if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
353 :
354 : /* now n and d are t_POLs in the same variable */
355 15841140 : v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
356 15841140 : if (!degpol(d))
357 : {
358 12764074 : n = RgX_Rg_div(n,gel(d,2));
359 12764074 : return v? RgX_mulXn(n,v): n;
360 : }
361 :
362 : /* X does not divide gcd(n,d), deg(d) > 0 */
363 3077066 : if (!isinexact(n) && !isinexact(d))
364 : {
365 3076821 : y = RgX_divrem(n, d, &z);
366 3076821 : if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
367 248038 : z = RgX_gcd(d, z);
368 248038 : if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
369 : }
370 248283 : return fix_rfrac(gred_rfrac_simple(n,d), v);
371 : }
372 :
373 : /* x,y t_INT, return x/y in reduced form */
374 : GEN
375 113534026 : Qdivii(GEN x, GEN y)
376 : {
377 113534026 : pari_sp av = avma;
378 : GEN r, q;
379 :
380 113534026 : if (lgefint(y) == 3)
381 : {
382 98169330 : q = Qdiviu(x, y[2]);
383 98165078 : if (signe(y) > 0) return q;
384 10672894 : if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
385 10673121 : return q;
386 : }
387 15364696 : if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
388 15365615 : if (equali1(x))
389 : {
390 5177514 : if (!signe(y)) pari_err_INV("gdiv",y);
391 5177402 : retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
392 : }
393 10188111 : q = dvmdii(x,y,&r);
394 10187931 : if (r == gen_0) return q; /* gen_0 intended */
395 7229962 : r = gcdii(y, r);
396 7229871 : if (lgefint(r) == 3)
397 : {
398 6392637 : ulong t = r[2];
399 6392637 : set_avma(av);
400 6392651 : if (t == 1) q = mkfraccopy(x,y);
401 : else
402 : {
403 2529239 : q = cgetg(3,t_FRAC);
404 2529446 : gel(q,1) = diviuexact(x,t);
405 2529310 : gel(q,2) = diviuexact(y,t);
406 : }
407 : }
408 : else
409 : { /* rare: r and q left on stack for efficiency */
410 837234 : q = cgetg(3,t_FRAC);
411 837257 : gel(q,1) = diviiexact(x,r);
412 837275 : gel(q,2) = diviiexact(y,r);
413 : }
414 7229927 : normalize_frac(q); return q;
415 : }
416 :
417 : /* x t_INT, return x/y in reduced form */
418 : GEN
419 116780366 : Qdiviu(GEN x, ulong y)
420 : {
421 116780366 : pari_sp av = avma;
422 : ulong r, t;
423 : GEN q;
424 :
425 116780366 : if (y == 1) return icopy(x);
426 95750645 : if (!y) pari_err_INV("gdiv",gen_0);
427 95750645 : if (equali1(x)) retmkfrac(gen_1, utoipos(y));
428 90037541 : q = absdiviu_rem(x,y,&r);
429 90030335 : if (!r)
430 : {
431 49584317 : if (signe(x) < 0) togglesign(q);
432 49584626 : return q;
433 : }
434 40446018 : t = ugcd(y, r); set_avma(av);
435 40451137 : if (t == 1) retmkfrac(icopy(x), utoipos(y));
436 16256189 : retmkfrac(diviuexact(x,t), utoipos(y / t));
437 : }
438 :
439 : /* x t_INT, return x/y in reduced form */
440 : GEN
441 1434225 : Qdivis(GEN x, long y)
442 : {
443 1434225 : pari_sp av = avma;
444 : ulong r, t;
445 : long s;
446 : GEN q;
447 :
448 1434225 : if (y > 0) return Qdiviu(x, y);
449 98745 : if (!y) pari_err_INV("gdiv",gen_0);
450 98745 : s = signe(x);
451 98745 : if (!s) return gen_0;
452 70231 : if (y < 0) { y = -y; s = -s; }
453 70231 : if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
454 69951 : if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
455 67886 : q = absdiviu_rem(x,y,&r);
456 67886 : if (!r)
457 : {
458 38458 : if (s < 0) togglesign(q);
459 38458 : return q;
460 : }
461 29428 : t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
462 29428 : if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
463 29428 : gel(q,1) = x; setsigne(x, s);
464 29428 : gel(q,2) = utoipos(y); return q;
465 : }
466 :
467 : /*******************************************************************/
468 : /* */
469 : /* CONJUGATION */
470 : /* */
471 : /*******************************************************************/
472 : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
473 : static GEN
474 18189 : quad_polmod_conj(GEN x, GEN y)
475 : {
476 : GEN z, u, v, a, b;
477 18189 : if (typ(x) != t_POL) return gcopy(x);
478 18189 : if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
479 18189 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
480 18189 : b = gel(y,3); v = gel(x,2);
481 18189 : z = cgetg(4, t_POL); z[1] = x[1];
482 18189 : gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
483 18189 : gel(z,3) = gneg(u); return z;
484 : }
485 : static GEN
486 18189 : quad_polmod_norm(GEN x, GEN y)
487 : {
488 : GEN z, u, v, a, b, c;
489 18189 : if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
490 18189 : a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
491 18189 : b = gel(y,3); v = gel(x,2);
492 18189 : c = gel(y,2);
493 18189 : z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
494 18189 : if (!gequal1(a)) z = gdiv(z, a);
495 18189 : return gadd(z, gsqr(v));
496 : }
497 :
498 : GEN
499 17711111 : conj_i(GEN x)
500 : {
501 : long lx, i;
502 : GEN y;
503 :
504 17711111 : switch(typ(x))
505 : {
506 5791189 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
507 5791189 : return x;
508 :
509 11759562 : case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
510 :
511 924 : case t_QUAD:
512 924 : y = cgetg(4,t_QUAD);
513 924 : gel(y,1) = gel(x,1);
514 924 : gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
515 924 : : gadd(gel(x,2), gel(x,3));
516 924 : gel(y,3) = gneg(gel(x,3)); return y;
517 :
518 9135 : case t_POL: case t_SER:
519 9135 : y = cgetg_copy(x, &lx); y[1] = x[1];
520 30597 : for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
521 9135 : return y;
522 :
523 150321 : case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
524 150321 : y = cgetg_copy(x, &lx);
525 542804 : for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
526 150324 : return y;
527 :
528 0 : case t_POLMOD:
529 : {
530 0 : GEN X = gel(x,1);
531 0 : long d = degpol(X);
532 0 : if (d < 2) return x;
533 0 : if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
534 : }
535 : }
536 0 : pari_err_TYPE("gconj",x);
537 : return NULL; /* LCOV_EXCL_LINE */
538 : }
539 : GEN
540 99240 : gconj(GEN x)
541 99240 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
542 :
543 : GEN
544 84 : conjvec(GEN x,long prec)
545 : {
546 : long lx, s, i;
547 : GEN z;
548 :
549 84 : switch(typ(x))
550 : {
551 0 : case t_INT: case t_INTMOD: case t_FRAC:
552 0 : return mkcolcopy(x);
553 :
554 0 : case t_COMPLEX: case t_QUAD:
555 0 : z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
556 :
557 28 : case t_FFELT:
558 28 : return FF_conjvec(x);
559 :
560 0 : case t_VEC: case t_COL:
561 0 : lx = lg(x); z = cgetg(lx,t_MAT);
562 0 : if (lx == 1) return z;
563 0 : gel(z,1) = conjvec(gel(x,1),prec);
564 0 : s = lgcols(z);
565 0 : for (i=2; i<lx; i++)
566 : {
567 0 : gel(z,i) = conjvec(gel(x,i),prec);
568 0 : if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
569 : }
570 0 : break;
571 :
572 56 : case t_POLMOD: {
573 56 : GEN T = gel(x,1), r;
574 : pari_sp av;
575 :
576 56 : lx = lg(T);
577 56 : if (lx <= 3) return cgetg(1,t_COL);
578 56 : x = gel(x,2);
579 238 : for (i=2; i<lx; i++)
580 : {
581 189 : GEN c = gel(T,i);
582 189 : switch(typ(c)) {
583 7 : case t_INTMOD: {
584 7 : GEN p = gel(c,1);
585 : pari_sp av;
586 7 : if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
587 7 : av = avma;
588 7 : T = RgX_to_FpX(T,p);
589 7 : x = RgX_to_FpX(x, p);
590 7 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
591 7 : z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
592 7 : return gerepileupto(av, z);
593 : }
594 182 : case t_INT:
595 182 : case t_FRAC: break;
596 0 : default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
597 : }
598 : }
599 49 : if (typ(x) != t_POL)
600 : {
601 0 : if (!is_rational_t(typ(x)))
602 0 : pari_err_TYPE("conjvec [not a rational t_POL]",x);
603 0 : retconst_col(lx-3, gcopy(x));
604 : }
605 49 : RgX_check_QX(x,"conjvec");
606 49 : av = avma;
607 49 : if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
608 49 : r = cleanroots(T,prec);
609 49 : z = cgetg(lx-2,t_COL);
610 182 : for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
611 49 : return gerepileupto(av, z);
612 : }
613 :
614 0 : default:
615 0 : pari_err_TYPE("conjvec",x);
616 : return NULL; /* LCOV_EXCL_LINE */
617 : }
618 0 : return z;
619 : }
620 :
621 : /********************************************************************/
622 : /** **/
623 : /** ADDITION **/
624 : /** **/
625 : /********************************************************************/
626 : /* x, y compatible PADIC, op = add or sub */
627 : static GEN
628 10556463 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
629 : {
630 10556463 : pari_sp av = avma;
631 : long d,e,r,rx,ry;
632 10556463 : GEN u, z, mod, p = gel(x,2);
633 : int swap;
634 :
635 10556463 : (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
636 10556108 : e = valp(x);
637 10556108 : r = valp(y); d = r-e;
638 10556108 : if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
639 10556108 : rx = precp(x);
640 10556108 : ry = precp(y);
641 10556108 : if (d) /* v(x) < v(y) */
642 : {
643 7386991 : r = d+ry; z = powiu(p,d);
644 7387124 : if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
645 7387136 : z = mulii(z,gel(y,4));
646 7387067 : u = swap? op(z, gel(x,4)): op(gel(x,4), z);
647 : }
648 : else
649 : {
650 : long c;
651 3169117 : if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
652 3169117 : u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
653 3170436 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
654 : {
655 138635 : set_avma(av); return zeropadic(p, e+r);
656 : }
657 3031662 : if (c)
658 : {
659 1874935 : mod = diviiexact(mod, powiu(p,c));
660 1874933 : r -= c;
661 1874933 : e += c;
662 : }
663 : }
664 10418707 : u = modii(u, mod);
665 10417396 : set_avma(av); z = cgetg(5,t_PADIC);
666 10416545 : z[1] = evalprecp(r) | evalvalp(e);
667 10416456 : gel(z,2) = icopy(p);
668 10416082 : gel(z,3) = icopy(mod);
669 10415719 : gel(z,4) = icopy(u); return z;
670 : }
671 : /* Rg_to_Fp(t_FRAC) without GC */
672 : static GEN
673 23559 : Q_to_Fp(GEN x, GEN p)
674 23559 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
675 : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
676 : static GEN
677 3515236 : addQp(GEN x, GEN y)
678 : {
679 3515236 : pari_sp av = avma;
680 3515236 : long d, r, e, vy = valp(y), py = precp(y);
681 3515236 : GEN z, q, mod, u, p = gel(y,2);
682 :
683 3515236 : e = Q_pvalrem(x, p, &x);
684 3515224 : d = vy - e; r = d + py;
685 3515224 : if (r <= 0) { set_avma(av); return gcopy(y); }
686 3513264 : mod = gel(y,3);
687 3513264 : u = gel(y,4);
688 3513264 : (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
689 :
690 3513253 : if (d > 0)
691 : {
692 893534 : q = powiu(p,d);
693 893620 : mod = mulii(mod, q);
694 893617 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
695 893617 : u = addii(x, mulii(u, q));
696 : }
697 2619719 : else if (d < 0)
698 : {
699 399310 : q = powiu(p,-d);
700 399310 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
701 399310 : u = addii(u, mulii(x, q));
702 399310 : r = py; e = vy;
703 : }
704 : else
705 : {
706 : long c;
707 2220409 : if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
708 2220409 : u = addii(u, x);
709 2220415 : if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
710 : {
711 906 : set_avma(av); return zeropadic(p,e+r);
712 : }
713 2219508 : if (c)
714 : {
715 797977 : mod = diviiexact(mod, powiu(p,c));
716 797977 : r -= c;
717 797977 : e += c;
718 : }
719 : }
720 3512444 : u = modii(u, mod); set_avma(av);
721 3512386 : z = cgetg(5,t_PADIC);
722 3512326 : z[1] = evalprecp(r) | evalvalp(e);
723 3512324 : gel(z,2) = icopy(p);
724 3512305 : gel(z,3) = icopy(mod);
725 3512275 : gel(z,4) = icopy(u); return z;
726 : }
727 :
728 : /* Mod(x,X) + Mod(y,X) */
729 : #define addsub_polmod_same addsub_polmod_scal
730 : /* Mod(x,X) +/- Mod(y,Y) */
731 : static GEN
732 2142 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
733 : {
734 2142 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
735 2142 : GEN z = cgetg(3,t_POLMOD);
736 2142 : long vx = varn(X), vy = varn(Y);
737 2142 : if (vx==vy) {
738 : pari_sp av;
739 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
740 14 : warn_coercion(X,Y,gel(z,1));
741 14 : gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
742 : }
743 2128 : if (varncmp(vx, vy) < 0)
744 2128 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
745 : else
746 0 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
747 2128 : gel(z,2) = op(x, y); return z;
748 : }
749 : /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
750 : static GEN
751 12531559 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
752 12531559 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
753 :
754 : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
755 : static GEN
756 396635 : add_ser_scal(GEN y, GEN x)
757 : {
758 : long i, v, ly, vy;
759 : GEN z;
760 :
761 396635 : if (isrationalzero(x)) return gcopy(y);
762 371316 : ly = lg(y);
763 371316 : v = valser(y);
764 371316 : if (v < 3-ly) return gcopy(y);
765 : /* v + ly >= 3 */
766 371092 : if (v < 0)
767 : {
768 1169 : z = cgetg(ly,t_SER); z[1] = y[1];
769 3304 : for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
770 1169 : gel(z,i) = gadd(x,gel(y,i)); i++;
771 3276 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
772 1169 : return normalizeser(z);
773 : }
774 369923 : vy = varn(y);
775 369923 : if (v > 0)
776 : {
777 16436 : if (ser_isexactzero(y))
778 7770 : return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
779 8666 : y -= v; ly += v;
780 8666 : z = cgetg(ly,t_SER);
781 8666 : x = gcopy(x);
782 19264 : for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
783 : }
784 353487 : else if (ser_isexactzero(y)) return gcopy(y);
785 : else
786 : { /* v = 0, ly >= 3 */
787 353480 : z = cgetg(ly,t_SER);
788 353480 : x = gadd(x, gel(y,2));
789 353480 : i = 3;
790 : }
791 1588493 : for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
792 362146 : gel(z,2) = x;
793 362146 : z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
794 362146 : return gequal0(x)? normalizeser(z): z;
795 : }
796 : static long
797 7063558 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
798 : /* x,y t_SER in the same variable: x+y */
799 : static GEN
800 3532178 : ser_add(GEN x, GEN y)
801 : {
802 3532178 : long i, lx,ly, n = valser(y) - valser(x);
803 : GEN z;
804 3532178 : if (n < 0) { n = -n; swap(x,y); }
805 : /* valser(x) <= valser(y) */
806 3532178 : lx = _serprec(x);
807 3532178 : if (lx == 2) /* don't lose type information */
808 : {
809 798 : z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
810 798 : setvalser(z, valser(x)); return z;
811 : }
812 3531380 : ly = _serprec(y) + n; if (lx < ly) ly = lx;
813 3531380 : if (n)
814 : {
815 107138 : if (n+2 > lx) return gcopy(x);
816 106564 : z = cgetg(ly,t_SER);
817 816687 : for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
818 505063 : for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
819 : } else {
820 3424242 : z = cgetg(ly,t_SER);
821 17326973 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
822 : }
823 3530806 : z[1] = x[1]; return normalizeser(z);
824 : }
825 : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
826 : static GEN
827 8830980 : add_rfrac_scal(GEN y, GEN x)
828 : {
829 : pari_sp av;
830 : GEN n;
831 :
832 8830980 : if (isintzero(x)) return gcopy(y); /* frequent special case */
833 5157173 : av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
834 5157173 : return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
835 : }
836 :
837 : /* x "scalar", ty != t_MAT and nonscalar */
838 : static GEN
839 22277448 : add_scal(GEN y, GEN x, long ty)
840 : {
841 22277448 : switch(ty)
842 : {
843 17381608 : case t_POL: return RgX_Rg_add(y, x);
844 396607 : case t_SER: return add_ser_scal(y, x);
845 4468813 : case t_RFRAC: return add_rfrac_scal(y, x);
846 0 : case t_COL: return RgC_Rg_add(y, x);
847 30414 : case t_VEC:
848 30414 : if (isintzero(x)) return gcopy(y);
849 182 : break;
850 : }
851 188 : pari_err_TYPE2("+",x,y);
852 : return NULL; /* LCOV_EXCL_LINE */
853 : }
854 :
855 : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
856 : static GEN
857 13373753 : setfrac(GEN z, GEN a, GEN b)
858 : {
859 13373753 : gel(z,1) = icopy_avma(a, (pari_sp)z);
860 13373809 : gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
861 13373864 : set_avma((pari_sp)gel(z,2)); return z;
862 : }
863 : /* z <- a / (b*Q), (Q,a) = 1 */
864 : static GEN
865 12408763 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
866 : {
867 12408763 : GEN q = Qdivii(a, b); /* != 0 */
868 12408870 : if (typ(q) == t_INT)
869 : {
870 1562465 : gel(z,1) = gerepileuptoint((pari_sp)Q, q);
871 1562465 : gel(z,2) = Q; return z;
872 : }
873 10846405 : return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
874 : }
875 : static GEN
876 25820437 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
877 : {
878 25820437 : GEN x1 = gel(x,1), x2 = gel(x,2);
879 25820437 : GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
880 25820437 : int s = cmpii(x2, y2);
881 :
882 : /* common denominator: (x1 op y1) / x2 */
883 25820446 : if (!s)
884 : {
885 9340046 : pari_sp av = avma;
886 9340046 : return gerepileupto(av, Qdivii(op(x1, y1), x2));
887 : }
888 16480400 : z = cgetg(3, t_FRAC);
889 16480425 : if (s < 0)
890 : {
891 9364131 : Q = dvmdii(y2, x2, &r);
892 : /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
893 9364040 : if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
894 2706978 : delta = gcdii(x2,r);
895 : }
896 : else
897 : {
898 7116294 : Q = dvmdii(x2, y2, &r);
899 : /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
900 7116303 : if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
901 1364522 : delta = gcdii(y2,r);
902 : }
903 : /* delta = gcd(x2,y2) */
904 4071571 : if (equali1(delta))
905 : { /* numerator is nonzero */
906 1544072 : gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
907 1544073 : gel(z,2) = mulii(x2,y2); return z;
908 : }
909 2527491 : x2 = diviiexact(x2,delta);
910 2527497 : y2 = diviiexact(y2,delta);
911 2527497 : n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
912 2527496 : q = dvmdii(n, delta, &r);
913 2527498 : if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
914 2316600 : r = gcdii(delta, r);
915 2316600 : if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
916 2316602 : return setfrac(z, n, mulii(mulii(x2, y2), delta));
917 : }
918 :
919 : /* assume x2, y2 are t_POLs in the same variable */
920 : static GEN
921 3039481 : add_rfrac(GEN x, GEN y)
922 : {
923 3039481 : pari_sp av = avma;
924 3039481 : GEN x1 = gel(x,1), x2 = gel(x,2);
925 3039481 : GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
926 :
927 3039481 : delta = RgX_gcd(x2,y2);
928 3039481 : if (!degpol(delta))
929 : {
930 672 : n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
931 672 : d = RgX_mul(x2, y2);
932 672 : return gerepileupto(av, gred_rfrac_simple(n, d));
933 : }
934 3038809 : x2 = RgX_div(x2,delta);
935 3038809 : y2 = RgX_div(y2,delta);
936 3038809 : n = gadd(gmul(x1,y2), gmul(y1,x2));
937 3038809 : if (!signe(n))
938 : {
939 721425 : n = simplify_shallow(n);
940 721425 : if (isexactzero(n))
941 : {
942 721418 : if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
943 14 : return gerepileupto(av, scalarpol(n, varn(x2)));
944 : }
945 7 : return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
946 : }
947 2317384 : if (degpol(n) == 0)
948 1150596 : return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
949 1166788 : q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
950 1166788 : if (isexactzero(r))
951 : {
952 : GEN z;
953 228094 : d = RgX_mul(x2, y2);
954 : /* "constant" denominator ? */
955 228094 : z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
956 228094 : return gerepileupto(av, z);
957 : }
958 938694 : r = RgX_gcd(delta, r);
959 938694 : if (degpol(r))
960 : {
961 160726 : n = RgX_div(n, r);
962 160726 : d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
963 : }
964 : else
965 777968 : d = RgX_mul(gel(x,2), y2);
966 938694 : return gerepileupto(av, gred_rfrac_simple(n, d));
967 : }
968 :
969 : GEN
970 4685088288 : gadd(GEN x, GEN y)
971 : {
972 4685088288 : long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
973 : pari_sp av, tetpil;
974 : GEN z, p1;
975 :
976 4685088288 : if (tx == ty) switch(tx) /* shortcut to generic case */
977 : {
978 1920646708 : case t_INT: return addii(x,y);
979 1508231311 : case t_REAL: return addrr(x,y);
980 1521938 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
981 1521938 : z = cgetg(3,t_INTMOD);
982 1521938 : if (X==Y || equalii(X,Y))
983 1521924 : return add_intmod_same(z, X, gel(x,2), gel(y,2));
984 14 : gel(z,1) = gcdii(X,Y);
985 14 : warn_coercion(X,Y,gel(z,1));
986 14 : av = avma; p1 = addii(gel(x,2),gel(y,2));
987 14 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
988 : }
989 20744706 : case t_FRAC: return addsub_frac(x,y,addii);
990 303492651 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
991 303698026 : gel(z,2) = gadd(gel(x,2),gel(y,2));
992 303579818 : if (isintzero(gel(z,2)))
993 : {
994 363265 : set_avma((pari_sp)(z+3));
995 363265 : return gadd(gel(x,1),gel(y,1));
996 : }
997 303088824 : gel(z,1) = gadd(gel(x,1),gel(y,1));
998 303241784 : return z;
999 2773407 : case t_PADIC:
1000 2773407 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
1001 2773229 : return addsub_pp(x,y, addii);
1002 672 : case t_QUAD: z = cgetg(4,t_QUAD);
1003 672 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1004 672 : gel(z,1) = ZX_copy(gel(x,1));
1005 672 : gel(z,2) = gadd(gel(x,2),gel(y,2));
1006 672 : gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
1007 9859682 : case t_POLMOD:
1008 9859682 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1009 9857575 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
1010 2107 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
1011 14006036 : case t_FFELT: return FF_add(x,y);
1012 66982910 : case t_POL:
1013 66982910 : vx = varn(x);
1014 66982910 : vy = varn(y);
1015 66982910 : if (vx != vy) {
1016 872095 : if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
1017 65742 : else return RgX_Rg_add(y, x);
1018 : }
1019 66110815 : return RgX_add(x, y);
1020 3526928 : case t_SER:
1021 3526928 : vx = varn(x);
1022 3526928 : vy = varn(y);
1023 3526928 : if (vx != vy) {
1024 28 : if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1025 21 : else return add_ser_scal(y, x);
1026 : }
1027 3526900 : return ser_add(x, y);
1028 4331382 : case t_RFRAC:
1029 4331382 : vx = varn(gel(x,2));
1030 4331382 : vy = varn(gel(y,2));
1031 4331382 : if (vx != vy) {
1032 1291901 : if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1033 538397 : else return add_rfrac_scal(y, x);
1034 : }
1035 3039481 : return add_rfrac(x,y);
1036 407348 : case t_VEC:
1037 407348 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1038 407348 : return RgV_add(x,y);
1039 1098491 : case t_COL:
1040 1098491 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1041 1098491 : return RgC_add(x,y);
1042 671307 : case t_MAT:
1043 671307 : lx = lg(x);
1044 671307 : if (lg(y) != lx) pari_err_OP("+",x,y);
1045 671307 : if (lx == 1) return cgetg(1, t_MAT);
1046 671307 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1047 671300 : return RgM_add(x,y);
1048 :
1049 0 : default: pari_err_TYPE2("+",x,y);
1050 : }
1051 : /* tx != ty */
1052 831127515 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1053 :
1054 831127515 : if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1055 : {
1056 746911990 : case t_INT:
1057 : switch(ty)
1058 : {
1059 479484235 : case t_REAL: return addir(x,y);
1060 2083559 : case t_INTMOD:
1061 2083559 : z = cgetg(3, t_INTMOD);
1062 2083559 : return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1063 44301128 : case t_FRAC: z = cgetg(3,t_FRAC);
1064 44301120 : gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1065 44301092 : gel(z,2) = icopy(gel(y,2)); return z;
1066 215081993 : case t_COMPLEX: return addRc(x, y);
1067 4643325 : case t_PADIC:
1068 4643325 : if (!signe(x)) return gcopy(y);
1069 3489535 : return addQp(x,y);
1070 1113 : case t_QUAD: return addRq(x, y);
1071 1356791 : case t_FFELT: return FF_Z_add(y,x);
1072 : }
1073 :
1074 : case t_REAL:
1075 49351638 : switch(ty)
1076 : {
1077 13791648 : case t_FRAC:
1078 13791648 : if (!signe(gel(y,1))) return rcopy(x);
1079 13791648 : if (!signe(x))
1080 : {
1081 11666 : lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1082 11666 : return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1083 : }
1084 13779982 : av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
1085 13778611 : return gerepile(av,tetpil,divri(z,gel(y,2)));
1086 35559920 : case t_COMPLEX: return addRc(x, y);
1087 70 : case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
1088 :
1089 0 : default: pari_err_TYPE2("+",x,y);
1090 : }
1091 :
1092 17640 : case t_INTMOD:
1093 : switch(ty)
1094 : {
1095 17507 : case t_FRAC: { GEN X = gel(x,1);
1096 17507 : z = cgetg(3, t_INTMOD);
1097 17507 : p1 = Fp_div(gel(y,1), gel(y,2), X);
1098 17507 : return add_intmod_same(z, X, p1, gel(x,2));
1099 : }
1100 14 : case t_FFELT:
1101 14 : if (!equalii(gel(x,1),FF_p_i(y)))
1102 0 : pari_err_OP("+",x,y);
1103 14 : return FF_Z_add(y,gel(x,2));
1104 91 : case t_COMPLEX: return addRc(x, y);
1105 0 : case t_PADIC: { GEN X = gel(x,1);
1106 0 : z = cgetg(3, t_INTMOD);
1107 0 : return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1108 : }
1109 28 : case t_QUAD: return addRq(x, y);
1110 : }
1111 :
1112 : case t_FRAC:
1113 : switch (ty)
1114 : {
1115 10527273 : case t_COMPLEX: return addRc(x, y);
1116 25694 : case t_PADIC:
1117 25694 : if (!signe(gel(x,1))) return gcopy(y);
1118 25694 : return addQp(x,y);
1119 133 : case t_QUAD: return addRq(x, y);
1120 1337 : case t_FFELT: return FF_Q_add(y, x);
1121 : }
1122 :
1123 : case t_FFELT:
1124 0 : pari_err_TYPE2("+",x,y);
1125 :
1126 35 : case t_COMPLEX:
1127 : switch(ty)
1128 : {
1129 28 : case t_PADIC:
1130 28 : return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
1131 7 : case t_QUAD:
1132 7 : lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1133 7 : return gequal0(y)? gcopy(x): addqf(y, x, lx);
1134 : }
1135 :
1136 : case t_PADIC: /* ty == t_QUAD */
1137 0 : return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
1138 : }
1139 : /* tx < ty, !is_const_t(y) */
1140 27448478 : switch(ty)
1141 : {
1142 8098 : case t_MAT:
1143 8098 : if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1144 8098 : if (isrationalzero(x)) return gcopy(y);
1145 8021 : return RgM_Rg_add(y, x);
1146 206828 : case t_COL:
1147 206828 : if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1148 206828 : return RgC_Rg_add(y, x);
1149 1839150 : case t_POLMOD: /* is_const_t(tx) in this case */
1150 1839150 : return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1151 : }
1152 25394402 : if (is_scalar_t(tx)) {
1153 22295637 : if (tx == t_POLMOD)
1154 : {
1155 107323 : vx = varn(gel(x,1));
1156 107323 : vy = gvar(y);
1157 107323 : if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1158 : else
1159 81486 : if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1160 40880 : return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1161 : }
1162 22188314 : return add_scal(y, x, ty);
1163 : }
1164 : /* x and y are not scalars, ty != t_MAT */
1165 3098758 : vx = gvar(x);
1166 3098758 : vy = gvar(y);
1167 3098758 : if (vx != vy) { /* x or y is treated as a scalar */
1168 22703 : if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1169 32298 : return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1170 32298 : : add_scal(y, x, ty);
1171 : }
1172 : /* vx = vy */
1173 3076055 : switch(tx)
1174 : {
1175 3075565 : case t_POL:
1176 : switch (ty)
1177 : {
1178 5299 : case t_SER:
1179 5299 : if (lg(x) == 2) return gcopy(y);
1180 5285 : i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1181 5285 : i = lg(y) + valser(y) - i;
1182 5285 : if (i < 3) return gcopy(y);
1183 5278 : p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1184 5278 : settyp(p1, t_VECSMALL); /* p1 left on stack */
1185 5278 : return y;
1186 :
1187 3070266 : case t_RFRAC: return add_rfrac_scal(y, x);
1188 : }
1189 0 : break;
1190 :
1191 490 : case t_SER:
1192 490 : if (ty == t_RFRAC)
1193 : {
1194 : long vn, vd;
1195 490 : av = avma;
1196 490 : vn = gval(gel(y,1), vy);
1197 490 : vd = RgX_valrem_inexact(gel(y,2), NULL);
1198 490 : if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
1199 :
1200 490 : l = lg(x) + valser(x) - (vn - vd);
1201 490 : if (l < 3) { set_avma(av); return gcopy(x); }
1202 490 : return gerepileupto(av, gadd(x, rfrac_to_ser_i(y, l)));
1203 : }
1204 0 : break;
1205 : }
1206 0 : pari_err_TYPE2("+",x,y);
1207 : return NULL; /* LCOV_EXCL_LINE */
1208 : }
1209 :
1210 : GEN
1211 287418056 : gaddsg(long x, GEN y)
1212 : {
1213 287418056 : long ty = typ(y);
1214 : GEN z;
1215 :
1216 287418056 : switch(ty)
1217 : {
1218 134410208 : case t_INT: return addsi(x,y);
1219 128770786 : case t_REAL: return addsr(x,y);
1220 11871 : case t_INTMOD:
1221 11871 : z = cgetg(3, t_INTMOD);
1222 11871 : return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1223 15247328 : case t_FRAC: z = cgetg(3,t_FRAC);
1224 15247328 : gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1225 15247328 : gel(z,2) = icopy(gel(y,2)); return z;
1226 8235576 : case t_COMPLEX:
1227 8235576 : z = cgetg(3, t_COMPLEX);
1228 8235576 : gel(z,1) = gaddsg(x, gel(y,1));
1229 8235576 : gel(z,2) = gcopy(gel(y,2)); return z;
1230 :
1231 742287 : default: return gadd(stoi(x), y);
1232 : }
1233 : }
1234 :
1235 : GEN
1236 3145598 : gsubsg(long x, GEN y)
1237 : {
1238 : GEN z, a, b;
1239 : pari_sp av;
1240 :
1241 3145598 : switch(typ(y))
1242 : {
1243 274876 : case t_INT: return subsi(x,y);
1244 1212965 : case t_REAL: return subsr(x,y);
1245 56 : case t_INTMOD:
1246 56 : z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1247 56 : return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1248 729169 : case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1249 729169 : gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
1250 729169 : gel(z,2) = icopy(gel(y,2)); return z;
1251 891663 : case t_COMPLEX:
1252 891663 : z = cgetg(3, t_COMPLEX);
1253 891663 : gel(z,1) = gsubsg(x, gel(y,1));
1254 891663 : gel(z,2) = gneg(gel(y,2)); return z;
1255 : }
1256 36869 : av = avma;
1257 36869 : return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
1258 : }
1259 :
1260 : /********************************************************************/
1261 : /** **/
1262 : /** SUBTRACTION **/
1263 : /** **/
1264 : /********************************************************************/
1265 :
1266 : GEN
1267 3652355822 : gsub(GEN x, GEN y)
1268 : {
1269 3652355822 : long tx = typ(x), ty = typ(y);
1270 : pari_sp av;
1271 : GEN z;
1272 3652355822 : if (tx == ty) switch(tx) /* shortcut to generic case */
1273 : {
1274 2911073182 : case t_INT: return subii(x,y);
1275 547075928 : case t_REAL: return subrr(x,y);
1276 1158483 : case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1277 1158483 : z = cgetg(3,t_INTMOD);
1278 1158483 : if (X==Y || equalii(X,Y))
1279 1158469 : return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1280 14 : gel(z,1) = gcdii(X,Y);
1281 14 : warn_coercion(X,Y,gel(z,1));
1282 14 : av = avma; p1 = subii(gel(x,2),gel(y,2));
1283 14 : gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
1284 : }
1285 5075761 : case t_FRAC: return addsub_frac(x,y, subii);
1286 110897159 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1287 110930498 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1288 110825650 : if (isintzero(gel(z,2)))
1289 : {
1290 21300 : set_avma((pari_sp)(z+3));
1291 21300 : return gsub(gel(x,1),gel(y,1));
1292 : }
1293 110808592 : gel(z,1) = gsub(gel(x,1),gel(y,1));
1294 110793967 : return z;
1295 7783299 : case t_PADIC:
1296 7783299 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
1297 7783299 : return addsub_pp(x,y, subii);
1298 168 : case t_QUAD: z = cgetg(4,t_QUAD);
1299 168 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1300 168 : gel(z,1) = ZX_copy(gel(x,1));
1301 168 : gel(z,2) = gsub(gel(x,2),gel(y,2));
1302 168 : gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1303 793989 : case t_POLMOD:
1304 793989 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1305 793954 : return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1306 35 : return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1307 203061 : case t_FFELT: return FF_sub(x,y);
1308 10248955 : case t_POL: {
1309 10248955 : long vx = varn(x);
1310 10248955 : long vy = varn(y);
1311 10248955 : if (vx != vy) {
1312 30263 : if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1313 5621 : else return Rg_RgX_sub(x, y);
1314 : }
1315 10218692 : return RgX_sub(x, y);
1316 : }
1317 297512 : case t_VEC:
1318 297512 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1319 297512 : return RgV_sub(x,y);
1320 7013101 : case t_COL:
1321 7013101 : if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1322 7013101 : return RgC_sub(x,y);
1323 235688 : case t_MAT: {
1324 235688 : long lx = lg(x);
1325 235688 : if (lg(y) != lx) pari_err_OP("+",x,y);
1326 235689 : if (lx == 1) return cgetg(1, t_MAT);
1327 156451 : if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1328 156451 : return RgM_sub(x,y);
1329 : }
1330 2458152 : case t_RFRAC: case t_SER: break;
1331 :
1332 0 : default: pari_err_TYPE2("+",x,y);
1333 : }
1334 49676871 : av = avma;
1335 52135023 : return gerepileupto(av, gadd(x,gneg_i(y)));
1336 : }
1337 :
1338 : /********************************************************************/
1339 : /** **/
1340 : /** MULTIPLICATION **/
1341 : /** **/
1342 : /********************************************************************/
1343 : static GEN
1344 306796 : mul_ser_scal(GEN y, GEN x)
1345 : {
1346 : long l, i;
1347 : GEN z;
1348 306796 : if (isexactzero(x)) return gmul(Rg_get_0(y), x);
1349 303548 : if (ser_isexactzero(y))
1350 : {
1351 518 : z = scalarser(lg(y) == 2? Rg_get_0(x): gmul(gel(y,2), x), varn(y), 1);
1352 518 : setvalser(z, valser(y)); return z;
1353 : }
1354 303030 : z = cgetg_copy(y, &l); z[1] = y[1];
1355 4182976 : for (i = 2; i < l; i++) gel(z,i) = gmul(gel(y,i), x);
1356 303030 : return normalizeser(z);
1357 : }
1358 : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1359 : * [n/d a valid RFRAC] */
1360 : static GEN
1361 10456316 : mul_rfrac_scal(GEN n, GEN d, GEN x)
1362 : {
1363 10456316 : pari_sp av = avma;
1364 : GEN z;
1365 :
1366 10456316 : switch(typ(x))
1367 : {
1368 21 : case t_PADIC:
1369 21 : n = gmul(n, x);
1370 21 : d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
1371 21 : return gerepileupto(av, gdiv(n,d));
1372 :
1373 504 : case t_INTMOD: case t_POLMOD:
1374 504 : n = gmul(n, x);
1375 504 : d = gmul(d, gmodulo(gen_1, gel(x,1)));
1376 504 : return gerepileupto(av, gdiv(n,d));
1377 : }
1378 10455791 : z = gred_rfrac2(x, d);
1379 10455791 : n = simplify_shallow(n);
1380 10455791 : if (typ(z) == t_RFRAC)
1381 : {
1382 7936847 : n = gmul(gel(z,1), n);
1383 7936847 : d = gel(z,2);
1384 7936847 : if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1385 0 : z = RgX_Rg_div(n, d);
1386 : else
1387 7936847 : z = gred_rfrac_simple(n, d);
1388 : }
1389 : else
1390 2518944 : z = gmul(z, n);
1391 10455791 : return gerepileupto(av, z);
1392 : }
1393 : static GEN
1394 112565159 : mul_scal(GEN y, GEN x, long ty)
1395 : {
1396 112565159 : switch(ty)
1397 : {
1398 103679755 : case t_POL:
1399 103679755 : if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1400 100522307 : return RgX_Rg_mul(y, x);
1401 298683 : case t_SER: return mul_ser_scal(y, x);
1402 8586739 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1403 14 : case t_QFB:
1404 14 : if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1405 : }
1406 12 : pari_err_TYPE2("*",x,y);
1407 : return NULL; /* LCOV_EXCL_LINE */
1408 : }
1409 :
1410 : static GEN
1411 160865 : mul_gen_rfrac(GEN X, GEN Y)
1412 : {
1413 160865 : GEN y1 = gel(Y,1), y2 = gel(Y,2);
1414 160865 : long vx = gvar(X), vy = varn(y2);
1415 166626 : return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1416 5761 : gred_rfrac_simple(gmul(y1,X), y2);
1417 : }
1418 : /* (x1/x2) * (y1/y2) */
1419 : static GEN
1420 7907114 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1421 : {
1422 : GEN z, X, Y;
1423 7907114 : pari_sp av = avma;
1424 :
1425 7907114 : X = gred_rfrac2(x1, y2);
1426 7907114 : Y = gred_rfrac2(y1, x2);
1427 7907114 : if (typ(X) == t_RFRAC)
1428 : {
1429 6628032 : if (typ(Y) == t_RFRAC) {
1430 6562456 : x1 = gel(X,1);
1431 6562456 : x2 = gel(X,2);
1432 6562456 : y1 = gel(Y,1);
1433 6562456 : y2 = gel(Y,2);
1434 6562456 : z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1435 : } else
1436 65576 : z = mul_gen_rfrac(Y, X);
1437 : }
1438 1279082 : else if (typ(Y) == t_RFRAC)
1439 95289 : z = mul_gen_rfrac(X, Y);
1440 : else
1441 1183793 : z = gmul(X, Y);
1442 7907114 : return gerepileupto(av, z);
1443 : }
1444 : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1445 : static GEN
1446 270821 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1447 : {
1448 270821 : pari_sp av = avma;
1449 270821 : GEN X = gred_rfrac2(x1, y2);
1450 270821 : if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1451 : {
1452 264255 : x2 = RgX_mul(gel(X,2), x2);
1453 264255 : x1 = gel(X,1);
1454 : }
1455 : else
1456 6566 : x1 = X;
1457 270821 : return gerepileupto(av, gred_rfrac_simple(x1, x2));
1458 : }
1459 :
1460 : /* Mod(y, Y) * x, assuming x scalar */
1461 : static GEN
1462 2578598 : mul_polmod_scal(GEN Y, GEN y, GEN x)
1463 2578598 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1464 :
1465 : /* cf mulqq */
1466 : static GEN
1467 5856548 : quad_polmod_mul(GEN P, GEN x, GEN y)
1468 : {
1469 5856548 : GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
1470 5856548 : pari_sp tetpil, av = avma;
1471 5856548 : T[1] = x[1];
1472 5856548 : p2 = gmul(gel(x,2), gel(y,2));
1473 5856548 : p3 = gmul(gel(x,3), gel(y,3));
1474 5856548 : p1 = gmul(gneg_i(c),p3);
1475 : /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
1476 5856548 : if (typ(b) == t_INT)
1477 : {
1478 5856527 : if (signe(b))
1479 : {
1480 4284755 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1481 4284755 : if (is_pm1(b))
1482 : {
1483 4284097 : if (signe(b) > 0) p3 = gneg(p3);
1484 : }
1485 : else
1486 658 : p3 = gmul(negi(b), p3);
1487 : }
1488 : else
1489 : {
1490 1571772 : p3 = gmul(gel(x,2),gel(y,3));
1491 1571772 : p4 = gmul(gel(x,3),gel(y,2));
1492 : }
1493 : }
1494 : else
1495 : {
1496 21 : p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1497 21 : p3 = gmul(gneg_i(b), p3);
1498 : }
1499 5856548 : tetpil = avma;
1500 5856548 : gel(T,2) = gadd(p2, p1);
1501 5856548 : gel(T,3) = gadd(p4, p3);
1502 5856548 : gerepilecoeffssp(av,tetpil,T+2,2);
1503 5856548 : return normalizepol_lg(T,4);
1504 : }
1505 : /* Mod(x,T) * Mod(y,T) */
1506 : static GEN
1507 7787069 : mul_polmod_same(GEN T, GEN x, GEN y)
1508 : {
1509 7787069 : GEN z = cgetg(3,t_POLMOD), a;
1510 7787069 : long v = varn(T), lx = lg(x), ly = lg(y);
1511 7787069 : gel(z,1) = RgX_copy(T);
1512 : /* x * y mod T optimised */
1513 7787069 : if (typ(x) != t_POL || varn(x) != v || lx <= 3
1514 7167315 : || typ(y) != t_POL || varn(y) != v || ly <= 3)
1515 1490477 : a = gmul(x, y);
1516 : else
1517 : {
1518 6296592 : if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1519 5851396 : a = quad_polmod_mul(T, x, y);
1520 : else
1521 445196 : a = RgXQ_mul(x, y, gel(z,1));
1522 : }
1523 7787069 : gel(z,2) = a; return z;
1524 : }
1525 : static GEN
1526 47060 : sqr_polmod(GEN T, GEN x)
1527 : {
1528 47060 : GEN a, z = cgetg(3,t_POLMOD);
1529 47060 : gel(z,1) = RgX_copy(T);
1530 47060 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1531 14240 : a = gsqr(x);
1532 : else
1533 : {
1534 32820 : pari_sp av = avma;
1535 32820 : a = RgXQ_sqr(x, gel(z,1));
1536 32820 : a = gerepileupto(av, a);
1537 : }
1538 47060 : gel(z,2) = a; return z;
1539 : }
1540 : /* Mod(x,X) * Mod(y,Y) */
1541 : static GEN
1542 2678207 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1543 : {
1544 2678207 : long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1545 2678207 : long vx = varn(X), vy = varn(Y);
1546 2678207 : GEN z = cgetg(3,t_POLMOD);
1547 :
1548 2678207 : if (vx==vy) {
1549 : pari_sp av;
1550 14 : gel(z,1) = RgX_gcd(X,Y); av = avma;
1551 14 : warn_coercion(X,Y,gel(z,1));
1552 14 : gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
1553 14 : return z;
1554 : }
1555 2678193 : if (varncmp(vx, vy) < 0)
1556 414316 : { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1557 : else
1558 2263877 : { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1559 2678193 : gel(z,2) = gmul(x, y); return z;
1560 : }
1561 :
1562 : #if 0 /* used by 3M only */
1563 : /* set z = x+y and return 1 if x,y have the same sign
1564 : * set z = x-y and return 0 otherwise */
1565 : static int
1566 : did_add(GEN x, GEN y, GEN *z)
1567 : {
1568 : long tx = typ(x), ty = typ(y);
1569 : if (tx == ty) switch(tx)
1570 : {
1571 : case t_INT: *z = addii(x,y); return 1;
1572 : case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1573 : case t_REAL:
1574 : if (signe(x) == -signe(y))
1575 : { *z = subrr(x,y); return 0; }
1576 : else
1577 : { *z = addrr(x,y); return 1; }
1578 : }
1579 : if (tx == t_REAL) switch(ty)
1580 : {
1581 : case t_INT:
1582 : if (signe(x) == -signe(y))
1583 : { *z = subri(x,y); return 0; }
1584 : else
1585 : { *z = addri(x,y); return 1; }
1586 : case t_FRAC:
1587 : if (signe(x) == -signe(gel(y,1)))
1588 : { *z = gsub(x,y); return 0; }
1589 : else
1590 : { *z = gadd(x,y); return 1; }
1591 : }
1592 : else if (ty == t_REAL) switch(tx)
1593 : {
1594 : case t_INT:
1595 : if (signe(x) == -signe(y))
1596 : { *z = subir(x,y); return 0; }
1597 : else
1598 : { *z = addir(x,y); return 1; }
1599 : case t_FRAC:
1600 : if (signe(gel(x,1)) == -signe(y))
1601 : { *z = gsub(x,y); return 0; }
1602 : else
1603 : { *z = gadd(x,y); return 1; }
1604 : }
1605 : *z = gadd(x,y); return 1;
1606 : }
1607 : #endif
1608 : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1609 : static GEN
1610 11081698 : mulcIR(GEN x, GEN y)
1611 : {
1612 11081698 : GEN z = cgetg(3,t_COMPLEX);
1613 11081899 : pari_sp av = avma;
1614 11081899 : gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
1615 11082143 : gel(z,2) = gmul(y, gel(x,1));
1616 11081808 : return z;
1617 :
1618 : }
1619 : /* x,y COMPLEX */
1620 : static GEN
1621 251585207 : mulcc(GEN x, GEN y)
1622 : {
1623 251585207 : GEN xr = gel(x,1), xi = gel(x,2);
1624 251585207 : GEN yr = gel(y,1), yi = gel(y,2);
1625 : GEN p1, p2, p3, p4, z;
1626 : pari_sp tetpil, av;
1627 :
1628 251585207 : if (isintzero(xr))
1629 : {
1630 15397896 : if (isintzero(yr)) {
1631 7173448 : av = avma;
1632 7173448 : return gerepileupto(av, gneg(gmul(xi,yi)));
1633 : }
1634 8224257 : return mulcIR(y, xi);
1635 : }
1636 236186508 : if (isintzero(yr)) return mulcIR(x, yi);
1637 :
1638 233344947 : z = cgetg(3,t_COMPLEX); av = avma;
1639 : #if 0
1640 : /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1641 : * e.g. xr + xi if exponents differ */
1642 : if (did_add(xr, xi, &p3))
1643 : {
1644 : if (did_add(yr, yi, &p4)) {
1645 : /* R = xr*yr - xi*yi
1646 : * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1647 : p1 = gmul(xr,yr);
1648 : p2 = gmul(xi,yi); p2 = gneg(p2);
1649 : p3 = gmul(p3, p4);
1650 : p4 = gsub(p2, p1);
1651 : } else {
1652 : /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1653 : * I = xr*yi + xi*yr */
1654 : p1 = gmul(p3,p4);
1655 : p3 = gmul(xr,yi);
1656 : p4 = gmul(xi,yr);
1657 : p2 = gsub(p3, p4);
1658 : }
1659 : } else {
1660 : if (did_add(yr, yi, &p4)) {
1661 : /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1662 : * I = xr*yi +xi*yr */
1663 : p1 = gmul(p3,p4);
1664 : p3 = gmul(xr,yi);
1665 : p4 = gmul(xi,yr);
1666 : p2 = gsub(p4, p3);
1667 : } else {
1668 : /* R = xr*yr - xi*yi
1669 : * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1670 : p3 = gneg( gmul(p3, p4) );
1671 : p1 = gmul(xr,yr);
1672 : p2 = gmul(xi,yi);
1673 : p4 = gadd(p1, p2);
1674 :
1675 : p2 = gneg(p2);
1676 : }
1677 : }
1678 : tetpil = avma;
1679 : gel(z,1) = gadd(p1,p2);
1680 : gel(z,2) = gadd(p3,p4);
1681 : #else
1682 233361927 : if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1683 : { /* 3M formula */
1684 558863 : p3 = addii(xr,xi);
1685 558863 : p4 = addii(yr,yi);
1686 558863 : p1 = mulii(xr,yr);
1687 558863 : p2 = mulii(xi,yi);
1688 558863 : p3 = mulii(p3,p4);
1689 558863 : p4 = addii(p2,p1);
1690 558863 : tetpil = avma;
1691 558863 : gel(z,1) = subii(p1,p2);
1692 558863 : gel(z,2) = subii(p3,p4);
1693 558863 : if (!signe(gel(z,2)))
1694 113232 : return gerepileuptoint((pari_sp)(z+3), gel(z,1));
1695 : }
1696 : else
1697 : { /* naive 4M formula: avoid all loss of accuracy */
1698 232803064 : p1 = gmul(xr,yr);
1699 232732656 : p2 = gmul(xi,yi);
1700 232714544 : p3 = gmul(xr,yi);
1701 232715960 : p4 = gmul(xi,yr);
1702 232720577 : tetpil = avma;
1703 232720577 : gel(z,1) = gsub(p1,p2);
1704 232606195 : gel(z,2) = gadd(p3,p4);
1705 232611353 : if (isintzero(gel(z,2)))
1706 : {
1707 49063 : cgiv(gel(z,2));
1708 49063 : return gerepileupto((pari_sp)(z+3), gel(z,1));
1709 : }
1710 : }
1711 : #endif
1712 :
1713 233012057 : gerepilecoeffssp(av,tetpil, z+1,2); return z;
1714 : }
1715 : /* x,y PADIC */
1716 : static GEN
1717 16878836 : mulpp(GEN x, GEN y) {
1718 16878836 : long l = valp(x) + valp(y);
1719 : pari_sp av;
1720 : GEN z, t;
1721 16878836 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
1722 16878802 : if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
1723 16662201 : if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
1724 :
1725 16325470 : t = (precp(x) > precp(y))? y: x;
1726 16325470 : z = cgetp(t); setvalp(z,l); av = avma;
1727 16325524 : affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
1728 16325399 : set_avma(av); return z;
1729 : }
1730 : /* x,y QUAD */
1731 : static GEN
1732 1106 : mulqq(GEN x, GEN y) {
1733 1106 : GEN z = cgetg(4,t_QUAD);
1734 1106 : GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
1735 : pari_sp av, tetpil;
1736 1106 : if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
1737 :
1738 1106 : gel(z,1) = ZX_copy(P); av = avma;
1739 1106 : p2 = gmul(gel(x,2),gel(y,2));
1740 1106 : p3 = gmul(gel(x,3),gel(y,3));
1741 1106 : p1 = gmul(gneg_i(c),p3);
1742 :
1743 1106 : if (signe(b))
1744 987 : p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
1745 : else
1746 : {
1747 119 : p3 = gmul(gel(x,2),gel(y,3));
1748 119 : p4 = gmul(gel(x,3),gel(y,2));
1749 : }
1750 1106 : tetpil = avma;
1751 1106 : gel(z,2) = gadd(p2,p1);
1752 1106 : gel(z,3) = gadd(p4,p3);
1753 1106 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
1754 : }
1755 :
1756 : GEN
1757 10422356 : mulcxI(GEN x)
1758 : {
1759 : GEN z;
1760 10422356 : switch(typ(x))
1761 : {
1762 1575595 : case t_INT: case t_REAL: case t_FRAC:
1763 1575595 : return mkcomplex(gen_0, x);
1764 8847233 : case t_COMPLEX:
1765 8847233 : if (isintzero(gel(x,1))) return gneg(gel(x,2));
1766 8846674 : z = cgetg(3,t_COMPLEX);
1767 8847184 : gel(z,1) = gneg(gel(x,2));
1768 8847784 : gel(z,2) = gel(x,1); return z;
1769 48 : default:
1770 48 : return gmul(gen_I(), x);
1771 : }
1772 : }
1773 : GEN
1774 198167 : mulcxmI(GEN x)
1775 : {
1776 : GEN z;
1777 198167 : switch(typ(x))
1778 : {
1779 84926 : case t_INT: case t_REAL: case t_FRAC:
1780 84926 : return mkcomplex(gen_0, gneg(x));
1781 113178 : case t_COMPLEX:
1782 113178 : if (isintzero(gel(x,1))) return gel(x,2);
1783 110567 : z = cgetg(3,t_COMPLEX);
1784 110567 : gel(z,1) = gel(x,2);
1785 110567 : gel(z,2) = gneg(gel(x,1)); return z;
1786 63 : default:
1787 63 : return gmul(mkcomplex(gen_0, gen_m1), x);
1788 : }
1789 : }
1790 : /* x * I^k */
1791 : GEN
1792 268250 : mulcxpowIs(GEN x, long k)
1793 : {
1794 268250 : switch (k & 3)
1795 : {
1796 50426 : case 1: return mulcxI(x);
1797 47646 : case 2: return gneg(x);
1798 52054 : case 3: return mulcxmI(x);
1799 : }
1800 118124 : return x;
1801 : }
1802 :
1803 : static GEN
1804 5626022 : init_ser(long l, long v, long e)
1805 : {
1806 5626022 : GEN z = cgetg(l, t_SER);
1807 5626022 : z[1] = evalvalser(e) | evalvarn(v) | evalsigne(1); return z;
1808 : }
1809 :
1810 : /* fill in coefficients of t_SER z from coeffs of t_POL y */
1811 : static GEN
1812 5614557 : fill_ser(GEN z, GEN y)
1813 : {
1814 5614557 : long i, lx = lg(z), ly = lg(y);
1815 5614557 : if (ly >= lx) {
1816 25192809 : for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
1817 : } else {
1818 332748 : for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1819 237123 : for ( ; i < lx; i++) gel(z,i) = gen_0;
1820 : }
1821 5614557 : return normalizeser(z);
1822 : }
1823 : /* assume typ(x) = t_VEC */
1824 : static int
1825 84 : is_ext_qfr(GEN x)
1826 84 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
1827 168 : && typ(gel(x,2)) == t_REAL; }
1828 :
1829 : GEN
1830 8147283669 : gmul(GEN x, GEN y)
1831 : {
1832 : long tx, ty, lx, ly, vx, vy, i, l;
1833 : pari_sp av, tetpil;
1834 : GEN z, p1;
1835 :
1836 8147283669 : if (x == y) return gsqr(x);
1837 7273860324 : tx = typ(x); ty = typ(y);
1838 7273860324 : if (tx == ty) switch(tx)
1839 : {
1840 3519225574 : case t_INT: return mulii(x,y);
1841 1694357165 : case t_REAL: return mulrr(x,y);
1842 1656504 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1843 1656504 : z = cgetg(3,t_INTMOD);
1844 1656504 : if (X==Y || equalii(X,Y))
1845 1656469 : return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1846 35 : gel(z,1) = gcdii(X,Y);
1847 35 : warn_coercion(X,Y,gel(z,1));
1848 35 : av = avma; p1 = mulii(gel(x,2),gel(y,2));
1849 35 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
1850 : }
1851 21408398 : case t_FRAC:
1852 : {
1853 21408398 : GEN x1 = gel(x,1), x2 = gel(x,2);
1854 21408398 : GEN y1 = gel(y,1), y2 = gel(y,2);
1855 21408398 : z=cgetg(3,t_FRAC);
1856 21408397 : p1 = gcdii(x1, y2);
1857 21408389 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1858 21408388 : p1 = gcdii(x2, y1);
1859 21408386 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1860 21408387 : tetpil = avma;
1861 21408387 : gel(z,2) = mulii(x2,y2);
1862 21408382 : gel(z,1) = mulii(x1,y1);
1863 21408389 : fix_frac_if_int_GC(z,tetpil); return z;
1864 : }
1865 243470805 : case t_COMPLEX: return mulcc(x, y);
1866 10371099 : case t_PADIC: return mulpp(x, y);
1867 882 : case t_QUAD: return mulqq(x, y);
1868 14525286 : case t_FFELT: return FF_mul(x, y);
1869 10217308 : case t_POLMOD:
1870 10217308 : if (RgX_equal_var(gel(x,1), gel(y,1)))
1871 7539101 : return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1872 2678207 : return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1873 45034712 : case t_POL:
1874 45034712 : vx = varn(x);
1875 45034712 : vy = varn(y);
1876 45034712 : if (vx != vy) {
1877 4868039 : if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1878 2082733 : else return RgX_Rg_mul(y, x);
1879 : }
1880 40166673 : return RgX_mul(x, y);
1881 :
1882 4347426 : case t_SER: {
1883 4347426 : vx = varn(x);
1884 4347426 : vy = varn(y);
1885 4347426 : if (vx != vy) {
1886 3675 : if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1887 3675 : return mul_ser_scal(y, x);
1888 : }
1889 4343751 : lx = minss(lg(x), lg(y));
1890 4343751 : if (lx == 2) return zeroser(vx, valser(x)+valser(y));
1891 4343737 : av = avma; z = init_ser(lx, vx, valser(x)+valser(y));
1892 4343737 : x = ser2pol_i(x, lx);
1893 4343737 : y = ser2pol_i(y, lx);
1894 4343737 : y = RgXn_mul(x, y, lx-2);
1895 4343737 : return gerepilecopy(av, fill_ser(z,y));
1896 : }
1897 42 : case t_VEC:
1898 42 : if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
1899 : /* fall through, handle extended t_QFB */
1900 47804 : case t_QFB: return qfbcomp(x,y);
1901 6720458 : case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1902 2748133 : case t_MAT: return RgM_mul(x, y);
1903 :
1904 1631 : case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1905 1631 : z = cgetg_copy(x, &l);
1906 1631 : if (l != lg(y)) break;
1907 17325 : for (i=1; i<l; i++)
1908 : {
1909 15694 : long yi = y[i];
1910 15694 : if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1911 15694 : z[i] = x[yi];
1912 : }
1913 1631 : return z;
1914 :
1915 0 : default:
1916 0 : pari_err_TYPE2("*",x,y);
1917 : }
1918 : /* tx != ty */
1919 1702092757 : if (is_const_t(ty) && is_const_t(tx)) {
1920 1577293050 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
1921 1577293050 : switch(tx) {
1922 1468958932 : case t_INT:
1923 : switch(ty)
1924 : {
1925 1021708467 : case t_REAL: return signe(x)? mulir(x,y): gen_0;
1926 1279983 : case t_INTMOD:
1927 1279983 : z = cgetg(3, t_INTMOD);
1928 1279982 : return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1929 63802227 : case t_FRAC:
1930 63802227 : if (!signe(x)) return gen_0;
1931 46401775 : z=cgetg(3,t_FRAC);
1932 46401989 : p1 = gcdii(x,gel(y,2));
1933 46400599 : if (equali1(p1))
1934 : {
1935 27867450 : set_avma((pari_sp)z);
1936 27867452 : gel(z,2) = icopy(gel(y,2));
1937 27867608 : gel(z,1) = mulii(gel(y,1), x);
1938 : }
1939 : else
1940 : {
1941 18533490 : x = diviiexact(x,p1); tetpil = avma;
1942 18532824 : gel(z,2) = diviiexact(gel(y,2), p1);
1943 18532593 : gel(z,1) = mulii(gel(y,1), x);
1944 18533095 : fix_frac_if_int_GC(z,tetpil);
1945 : }
1946 46401404 : return z;
1947 377105780 : case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1948 5683858 : case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1949 1701 : case t_QUAD: return mulRq(x,y);
1950 1622392 : case t_FFELT: return FF_Z_mul(y,x);
1951 : }
1952 :
1953 : case t_REAL:
1954 100266599 : switch(ty)
1955 : {
1956 30207416 : case t_FRAC: return mulrfrac(x, y);
1957 70059141 : case t_COMPLEX: return mulRc(x, y);
1958 21 : case t_QUAD: return mulqf(y, x, realprec(x));
1959 21 : default: pari_err_TYPE2("*",x,y);
1960 : }
1961 :
1962 8015 : case t_INTMOD:
1963 : switch(ty)
1964 : {
1965 7133 : case t_FRAC: { GEN X = gel(x,1);
1966 7133 : z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1967 7133 : return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1968 : }
1969 49 : case t_COMPLEX: return mulRc_direct(x,y);
1970 420 : case t_PADIC: { GEN X = gel(x,1);
1971 420 : z = cgetg(3, t_INTMOD);
1972 420 : return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1973 : }
1974 63 : case t_QUAD: return mulRq(x, y);
1975 350 : case t_FFELT:
1976 350 : if (!equalii(gel(x,1),FF_p_i(y)))
1977 0 : pari_err_OP("*",x,y);
1978 350 : return FF_Z_mul(y,gel(x,2));
1979 : }
1980 :
1981 : case t_FRAC:
1982 : switch(ty)
1983 : {
1984 5203685 : case t_COMPLEX: return mulRc(x, y);
1985 3335477 : case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
1986 343 : case t_QUAD: return mulRq(x, y);
1987 1995 : case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
1988 : }
1989 :
1990 : case t_FFELT:
1991 36 : pari_err_TYPE2("*",x,y);
1992 :
1993 21 : case t_COMPLEX:
1994 : switch(ty)
1995 : {
1996 14 : case t_PADIC:
1997 14 : return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
1998 7 : case t_QUAD:
1999 7 : lx = precision(x); if (!lx) pari_err_OP("*",x,y);
2000 7 : return mulqf(y, x, lx);
2001 : }
2002 :
2003 : case t_PADIC: /* ty == t_QUAD */
2004 28 : return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
2005 : }
2006 : }
2007 :
2008 127298562 : if (is_matvec_t(ty))
2009 : {
2010 11338558 : if (!is_matvec_t(tx))
2011 : {
2012 10384144 : if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
2013 10384145 : z = cgetg_copy(y, &ly);
2014 1572393935 : for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
2015 10383828 : return z;
2016 : }
2017 954417 : switch(tx)
2018 : {
2019 121730 : case t_VEC:
2020 : switch(ty) {
2021 15540 : case t_COL: return RgV_RgC_mul(x,y);
2022 106190 : case t_MAT: return RgV_RgM_mul(x,y);
2023 : }
2024 0 : break;
2025 1687 : case t_COL:
2026 : switch(ty) {
2027 1687 : case t_VEC: return RgC_RgV_mul(x,y);
2028 0 : case t_MAT: return RgC_RgM_mul(x,y);
2029 : }
2030 0 : break;
2031 831014 : case t_MAT:
2032 : switch(ty) {
2033 0 : case t_VEC: return RgM_RgV_mul(x,y);
2034 831014 : case t_COL: return RgM_RgC_mul(x,y);
2035 : }
2036 : }
2037 : }
2038 118949128 : if (is_matvec_t(tx))
2039 : {
2040 1836877 : if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2041 1837074 : z = cgetg_copy(x, &lx);
2042 7584189 : for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
2043 1836980 : return z;
2044 : }
2045 117112490 : if (tx > ty) { swap(x,y); lswap(tx,ty); }
2046 : /* tx < ty, !ismatvec(x and y) */
2047 :
2048 117112490 : if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2049 2531215 : return mul_polmod_scal(gel(y,1), gel(y,2), x);
2050 114581275 : if (is_scalar_t(tx)) {
2051 111285985 : if (tx == t_POLMOD) {
2052 3100814 : vx = varn(gel(x,1));
2053 3100814 : vy = gvar(y);
2054 3100814 : if (vx != vy) {
2055 2853406 : if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2056 47383 : return mul_polmod_scal(gel(x,1), gel(x,2), y);
2057 : }
2058 : /* error if ty == t_SER */
2059 247408 : av = avma; y = gmod(y, gel(x,1));
2060 247401 : return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2061 : }
2062 108185171 : return mul_scal(y, x, ty);
2063 : }
2064 :
2065 : /* x and y are not scalars, nor matvec */
2066 3295134 : vx = gvar(x);
2067 3295129 : vy = gvar(y);
2068 3295129 : if (vx != vy) /* x or y is treated as a scalar */
2069 2780201 : return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2070 2780201 : : mul_scal(y, x, ty);
2071 : /* vx = vy */
2072 1876248 : switch(tx)
2073 : {
2074 1876220 : case t_POL:
2075 : switch (ty)
2076 : {
2077 6713 : case t_SER:
2078 : {
2079 : long v;
2080 6713 : av = avma; v = RgX_valrem(x, &x);
2081 6713 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2082 6699 : v += valser(y); ly = lg(y);
2083 6699 : if (ly == 2) { set_avma(av); return zeroser(vx, v); }
2084 6468 : if (degpol(x))
2085 : {
2086 2030 : z = init_ser(ly, vx, v);
2087 2030 : x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
2088 2030 : return gerepilecopy(av, fill_ser(z, x));
2089 : }
2090 : /* take advantage of x = c*t^v */
2091 4438 : set_avma(av); y = mul_ser_scal(y, gel(x,2));
2092 4438 : setvalser(y, v); return y;
2093 : }
2094 :
2095 1869507 : case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2096 : }
2097 0 : break;
2098 :
2099 28 : case t_SER:
2100 : switch (ty)
2101 : {
2102 28 : case t_RFRAC:
2103 28 : av = avma;
2104 28 : return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2105 : }
2106 0 : break;
2107 : }
2108 0 : pari_err_TYPE2("*",x,y);
2109 : return NULL; /* LCOV_EXCL_LINE */
2110 : }
2111 :
2112 : /* return a nonnormalized result */
2113 : GEN
2114 104212 : sqr_ser_part(GEN x, long l1, long l2)
2115 : {
2116 : long i, j, l;
2117 : pari_sp av;
2118 : GEN Z, z, p1, p2;
2119 : long mi;
2120 104212 : if (l2 < l1) return zeroser(varn(x), 2*valser(x));
2121 104198 : p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2122 104198 : Z = cgetg(l2-l1+3, t_SER);
2123 104198 : Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
2124 104198 : z = Z + 2-l1;
2125 104198 : x += 2; mi = 0;
2126 415856 : for (i=0; i<l1; i++)
2127 : {
2128 311658 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2129 : }
2130 :
2131 685212 : for (i=l1; i<=l2; i++)
2132 : {
2133 581014 : p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2134 581014 : p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2135 1217879 : for (j=i-mi; j<=minss(l,mi); j++)
2136 636865 : if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2137 581014 : p1 = gshift(p1,1);
2138 581014 : if ((i&1) == 0 && p2[i>>1])
2139 81069 : p1 = gadd(p1, gsqr(gel(x,i>>1)));
2140 581014 : gel(z,i) = gerepileupto(av,p1);
2141 : }
2142 104198 : return Z;
2143 : }
2144 :
2145 : GEN
2146 1074649937 : gsqr(GEN x)
2147 : {
2148 : long i, lx;
2149 : pari_sp av, tetpil;
2150 : GEN z, p1, p2, p3, p4;
2151 :
2152 1074649937 : switch(typ(x))
2153 : {
2154 883151795 : case t_INT: return sqri(x);
2155 177276383 : case t_REAL: return sqrr(x);
2156 142611 : case t_INTMOD: { GEN X = gel(x,1);
2157 142611 : z = cgetg(3,t_INTMOD);
2158 142611 : gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
2159 142613 : gel(z,1) = icopy(X); return z;
2160 : }
2161 3031814 : case t_FRAC: return sqrfrac(x);
2162 :
2163 6849824 : case t_COMPLEX:
2164 6849824 : if (isintzero(gel(x,1))) {
2165 210816 : av = avma;
2166 210816 : return gerepileupto(av, gneg(gsqr(gel(x,2))));
2167 : }
2168 6639001 : z = cgetg(3,t_COMPLEX); av = avma;
2169 6638989 : p1 = gadd(gel(x,1),gel(x,2));
2170 6638840 : p2 = gsub(gel(x,1), gel(x,2));
2171 6638788 : p3 = gmul(gel(x,1),gel(x,2));
2172 6638927 : tetpil = avma;
2173 6638927 : gel(z,1) = gmul(p1,p2);
2174 6638953 : gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
2175 :
2176 2597 : case t_PADIC:
2177 2597 : z = cgetg(5,t_PADIC);
2178 2597 : i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
2179 2597 : if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
2180 2597 : z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
2181 2597 : gel(z,2) = icopy(gel(x,2));
2182 2597 : gel(z,3) = shifti(gel(x,3), i); av = avma;
2183 2597 : gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
2184 2597 : return z;
2185 :
2186 70 : case t_QUAD: z = cgetg(4,t_QUAD);
2187 70 : p1 = gel(x,1);
2188 70 : gel(z,1) = ZX_copy(p1); av = avma;
2189 70 : p2 = gsqr(gel(x,2));
2190 70 : p3 = gsqr(gel(x,3));
2191 70 : p4 = gmul(gneg_i(gel(p1,2)),p3);
2192 :
2193 70 : if (gequal0(gel(p1,3)))
2194 : {
2195 7 : tetpil = avma;
2196 7 : gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
2197 7 : av = avma;
2198 7 : p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
2199 7 : gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
2200 : }
2201 :
2202 63 : p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
2203 63 : tetpil = avma;
2204 63 : gel(z,2) = gadd(p2,p4);
2205 63 : gel(z,3) = gadd(p1,p3);
2206 63 : gerepilecoeffssp(av,tetpil,z+2,2); return z;
2207 :
2208 47060 : case t_POLMOD:
2209 47060 : return sqr_polmod(gel(x,1), gel(x,2));
2210 :
2211 2940251 : case t_FFELT: return FF_sqr(x);
2212 :
2213 1107866 : case t_POL: return RgX_sqr(x);
2214 :
2215 26740 : case t_SER:
2216 26740 : lx = lg(x);
2217 26740 : if (ser_isexactzero(x)) {
2218 21 : GEN z = gcopy(x);
2219 21 : setvalser(z, 2*valser(x));
2220 21 : return z;
2221 : }
2222 26719 : if (lx < 40)
2223 26460 : return normalizeser( sqr_ser_part(x, 0, lx-3) );
2224 : else
2225 : {
2226 259 : pari_sp av = avma;
2227 259 : GEN z = init_ser(lx, varn(x), 2*valser(x));
2228 259 : x = ser2pol_i(x, lx);
2229 259 : x = RgXn_sqr(x, lx-2);
2230 259 : return gerepilecopy(av, fill_ser(z,x));
2231 : }
2232 :
2233 14 : case t_RFRAC: z = cgetg(3,t_RFRAC);
2234 14 : gel(z,1) = gsqr(gel(x,1));
2235 14 : gel(z,2) = gsqr(gel(x,2)); return z;
2236 :
2237 1099 : case t_MAT: return RgM_sqr(x);
2238 14 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
2239 : /* fall through handle extended t_QFB */
2240 81901 : case t_QFB: return qfbsqr(x);
2241 658 : case t_VECSMALL:
2242 658 : z = cgetg_copy(x, &lx);
2243 16289 : for (i=1; i<lx; i++)
2244 : {
2245 15631 : long xi = x[i];
2246 15631 : if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2247 15631 : z[i] = x[xi];
2248 : }
2249 658 : return z;
2250 : }
2251 7 : pari_err_TYPE2("*",x,x);
2252 : return NULL; /* LCOV_EXCL_LINE */
2253 : }
2254 :
2255 : /********************************************************************/
2256 : /** **/
2257 : /** DIVISION **/
2258 : /** **/
2259 : /********************************************************************/
2260 : static GEN
2261 30832 : div_rfrac_scal(GEN x, GEN y)
2262 : {
2263 30832 : pari_sp av = avma;
2264 30832 : GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2265 30832 : return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
2266 : }
2267 : static GEN
2268 37459 : div_scal_rfrac(GEN x, GEN y)
2269 : {
2270 37459 : GEN y1 = gel(y,1), y2 = gel(y,2);
2271 37459 : pari_sp av = avma;
2272 37459 : if (typ(y1) == t_POL && varn(y2) == varn(y1))
2273 : {
2274 7 : if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
2275 0 : y1 = gel(y1,2);
2276 : }
2277 37452 : return RgX_Rg_mul(y2, gdiv(x,y1));
2278 : }
2279 : static GEN
2280 1186656 : div_rfrac(GEN x, GEN y)
2281 1186656 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2282 :
2283 : /* x != 0 */
2284 : static GEN
2285 1321422 : div_ser_scal(GEN y, GEN x)
2286 : {
2287 : long i, l;
2288 : GEN z;
2289 1321422 : if (ser_isexactzero(y))
2290 : {
2291 28 : z = scalarser(lg(y) == 2? Rg_get_0(x): gdiv(gel(y,2), x), varn(y), 1);
2292 28 : setvalser(z, valser(y)); return z;
2293 : }
2294 1321394 : z = cgetg_copy(y, &l); z[1] = y[1];
2295 6108824 : for (i = 2; i < l; i++) gel(z,i) = gdiv(gel(y,i), x);
2296 1321394 : return normalizeser(z);
2297 : }
2298 : GEN
2299 7553 : ser_normalize(GEN x)
2300 : {
2301 7553 : long i, lx = lg(x);
2302 : GEN c, z;
2303 7553 : if (lx == 2) return x;
2304 7553 : c = gel(x,2); if (gequal1(c)) return x;
2305 7476 : z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2306 107835 : for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2307 7476 : return z;
2308 : }
2309 :
2310 : /* y != 0 */
2311 : static GEN
2312 5434393 : div_T_scal(GEN x, GEN y, long tx) {
2313 5434393 : switch(tx)
2314 : {
2315 4086013 : case t_POL: return RgX_Rg_div(x, y);
2316 1321415 : case t_SER: return div_ser_scal(x, y);
2317 26968 : case t_RFRAC: return div_rfrac_scal(x,y);
2318 : }
2319 0 : pari_err_TYPE2("/",x,y);
2320 : return NULL; /* LCOV_EXCL_LINE */
2321 : }
2322 :
2323 : static GEN
2324 9258912 : div_scal_pol(GEN x, GEN y) {
2325 9258912 : long ly = lg(y);
2326 : pari_sp av;
2327 9258912 : if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2328 9180637 : if (isrationalzero(x)) return zeropol(varn(y));
2329 7130767 : av = avma;
2330 7130767 : return gerepileupto(av, gred_rfrac_simple(x,y));
2331 : }
2332 : static GEN
2333 18382 : div_scal_ser(GEN x, GEN y)
2334 18382 : { pari_sp av = avma; return gerepileupto(av, gmul(x, ser_inv(y))); }
2335 : static GEN
2336 9267531 : div_scal_T(GEN x, GEN y, long ty) {
2337 9267531 : switch(ty)
2338 : {
2339 9212453 : case t_POL: return div_scal_pol(x, y);
2340 18375 : case t_SER: return div_scal_ser(x, y);
2341 36703 : case t_RFRAC: return div_scal_rfrac(x, y);
2342 : }
2343 0 : pari_err_TYPE2("/",x,y);
2344 : return NULL; /* LCOV_EXCL_LINE */
2345 : }
2346 :
2347 : /* assume tx = ty = t_SER, same variable vx */
2348 : static GEN
2349 761864 : div_ser(GEN x, GEN y, long vx)
2350 : {
2351 761864 : long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
2352 761864 : GEN y0 = y, z;
2353 761864 : pari_sp av = avma;
2354 :
2355 761864 : if (!signe(y)) pari_err_INV("div_ser", y);
2356 761857 : if (ser_isexactzero(x))
2357 : {
2358 59892 : if (lx == 2) return zeroser(vx, v);
2359 147 : z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
2360 147 : setvalser(z, v); return z;
2361 : }
2362 701965 : if (lx < ly) ly = lx;
2363 701965 : y = ser2pol_i_normalize(y, ly, &e);
2364 701965 : if (e)
2365 : {
2366 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2367 7 : ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
2368 : }
2369 701965 : z = init_ser(ly, vx, v);
2370 701965 : if (ly == 3)
2371 : {
2372 11465 : gel(z,2) = gdiv(gel(x,2), gel(y,2));
2373 11465 : if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
2374 11465 : return gerepileupto(av, z);
2375 : }
2376 690500 : x = ser2pol_i(x, ly);
2377 690500 : y = RgXn_div_i(x, y, ly-2);
2378 690500 : return gerepilecopy(av, fill_ser(z,y));
2379 : }
2380 : /* x,y compatible PADIC */
2381 : static GEN
2382 2824389 : divpp(GEN x, GEN y) {
2383 : pari_sp av;
2384 : long a, b;
2385 : GEN z, M;
2386 :
2387 2824389 : if (!signe(gel(y,4))) pari_err_INV("divpp",y);
2388 2824388 : if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
2389 2824073 : a = precp(x);
2390 2824073 : b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
2391 2824073 : z = cgetg(5, t_PADIC);
2392 2823761 : z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
2393 2823706 : gel(z,2) = icopy(gel(x,2));
2394 2823479 : gel(z,3) = icopy(M); av = avma;
2395 2823192 : gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
2396 2823746 : return z;
2397 : }
2398 : static GEN
2399 27831 : div_polmod_same(GEN T, GEN x, GEN y)
2400 : {
2401 27831 : long v = varn(T);
2402 27831 : GEN a, z = cgetg(3, t_POLMOD);
2403 27831 : gel(z,1) = RgX_copy(T);
2404 27831 : if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2405 10164 : a = gdiv(x, y);
2406 17667 : else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2407 119 : {
2408 119 : pari_sp av = avma;
2409 119 : a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
2410 : }
2411 17548 : else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2412 5152 : {
2413 5152 : pari_sp av = avma;
2414 5152 : a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2415 5152 : a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2416 5152 : a = gerepileupto(av, a);
2417 : }
2418 : else
2419 : {
2420 12396 : pari_sp av = avma;
2421 12396 : a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2422 12396 : a = gerepileupto(av, a);
2423 : }
2424 27831 : gel(z,2) = a; return z;
2425 : }
2426 : GEN
2427 417493564 : gdiv(GEN x, GEN y)
2428 : {
2429 417493564 : long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2430 : pari_sp av, tetpil;
2431 : GEN z, p1, p2;
2432 :
2433 417493564 : if (tx == ty) switch(tx)
2434 : {
2435 80332392 : case t_INT:
2436 80332392 : return Qdivii(x,y);
2437 :
2438 82194198 : case t_REAL: return divrr(x,y);
2439 25619 : case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2440 25619 : z = cgetg(3,t_INTMOD);
2441 25619 : if (X==Y || equalii(X,Y))
2442 25605 : return div_intmod_same(z, X, gel(x,2), gel(y,2));
2443 14 : gel(z,1) = gcdii(X,Y);
2444 14 : warn_coercion(X,Y,gel(z,1));
2445 14 : av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2446 14 : gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
2447 : }
2448 3753396 : case t_FRAC: {
2449 3753396 : GEN x1 = gel(x,1), x2 = gel(x,2);
2450 3753396 : GEN y1 = gel(y,1), y2 = gel(y,2);
2451 3753396 : z = cgetg(3, t_FRAC);
2452 3753396 : p1 = gcdii(x1, y1);
2453 3753394 : if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2454 3753395 : p1 = gcdii(x2, y2);
2455 3753394 : if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2456 3753392 : tetpil = avma;
2457 3753392 : gel(z,2) = mulii(x2,y1);
2458 3753392 : gel(z,1) = mulii(x1,y2);
2459 3753394 : normalize_frac(z);
2460 3753394 : fix_frac_if_int_GC(z,tetpil);
2461 3753396 : return z;
2462 : }
2463 8222089 : case t_COMPLEX:
2464 8222089 : if (isintzero(gel(y,1)))
2465 : {
2466 110495 : y = gel(y,2);
2467 110495 : if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2468 93415 : z = cgetg(3,t_COMPLEX);
2469 93415 : gel(z,1) = gdiv(gel(x,2), y);
2470 93415 : av = avma;
2471 93415 : gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
2472 93415 : return z;
2473 : }
2474 8111582 : av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
2475 8111672 : return gerepile(av, tetpil, gdiv(p2,p1));
2476 :
2477 775117 : case t_PADIC:
2478 775117 : if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
2479 775117 : return divpp(x, y);
2480 :
2481 238 : case t_QUAD:
2482 238 : if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2483 224 : av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
2484 224 : return gerepile(av, tetpil, gdiv(p2,p1));
2485 :
2486 236711 : case t_FFELT: return FF_div(x,y);
2487 :
2488 32206 : case t_POLMOD:
2489 32206 : if (RgX_equal_var(gel(x,1), gel(y,1)))
2490 27831 : z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2491 : else {
2492 4375 : av = avma;
2493 4375 : z = gerepileupto(av, gmul(x, ginv(y)));
2494 : }
2495 32206 : return z;
2496 :
2497 18461133 : case t_POL:
2498 18461133 : vx = varn(x);
2499 18461133 : vy = varn(y);
2500 18461133 : if (vx != vy) {
2501 102767 : if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2502 46459 : else return div_scal_pol(x, y);
2503 : }
2504 18358366 : if (!signe(y)) pari_err_INV("gdiv",y);
2505 18358366 : if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2506 18234870 : av = avma;
2507 18234870 : return gerepileupto(av, gred_rfrac2(x,y));
2508 :
2509 761885 : case t_SER:
2510 761885 : vx = varn(x);
2511 761885 : vy = varn(y);
2512 761885 : if (vx != vy) {
2513 21 : if (varncmp(vx, vy) < 0)
2514 : {
2515 14 : if (!signe(y)) pari_err_INV("gdiv",y);
2516 7 : return div_ser_scal(x, y);
2517 : }
2518 7 : return div_scal_ser(x, y);
2519 : }
2520 761864 : return div_ser(x, y, vx);
2521 1191136 : case t_RFRAC:
2522 1191136 : vx = varn(gel(x,2));
2523 1191136 : vy = varn(gel(y,2));
2524 1191136 : if (vx != vy) {
2525 4480 : if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2526 756 : else return div_scal_rfrac(x, y);
2527 : }
2528 1186656 : return div_rfrac(x,y);
2529 :
2530 21 : case t_VEC: /* handle extended t_QFB */
2531 21 : case t_QFB: av = avma; return gerepileupto(av, qfbcomp(x, ginv(y)));
2532 :
2533 14 : case t_MAT:
2534 14 : av = avma; p1 = RgM_inv(y);
2535 14 : if (!p1) pari_err_INV("gdiv",y);
2536 14 : return gerepileupto(av, RgM_mul(x, p1));
2537 :
2538 0 : default: pari_err_TYPE2("/",x,y);
2539 : }
2540 :
2541 221512581 : if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2542 : {
2543 3821086 : long s = signe(x);
2544 3821086 : if (!s) {
2545 479320 : if (gequal0(y)) pari_err_INV("gdiv",y);
2546 479320 : switch (ty)
2547 : {
2548 476156 : default: return gen_0;
2549 28 : case t_INTMOD:
2550 28 : z = cgetg(3,t_INTMOD);
2551 28 : gel(z,1) = icopy(gel(y,1));
2552 28 : gel(z,2) = gen_0; return z;
2553 3136 : case t_FFELT: return FF_zero(y);
2554 : }
2555 : }
2556 3341766 : if (is_pm1(x)) {
2557 1103399 : if (s > 0) return ginv(y);
2558 190340 : av = avma; return gerepileupto(av, ginv(gneg(y)));
2559 : }
2560 2238368 : switch(ty)
2561 : {
2562 782918 : case t_REAL: return divir(x,y);
2563 21 : case t_INTMOD:
2564 21 : z = cgetg(3, t_INTMOD);
2565 21 : return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2566 783374 : case t_FRAC:
2567 783374 : z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2568 783374 : if (equali1(p1))
2569 : {
2570 251166 : set_avma((pari_sp)z);
2571 251166 : gel(z,2) = icopy(gel(y,1));
2572 251166 : gel(z,1) = mulii(gel(y,2), x);
2573 251166 : normalize_frac(z);
2574 251166 : fix_frac_if_int(z);
2575 : }
2576 : else
2577 : {
2578 532208 : x = diviiexact(x,p1); tetpil = avma;
2579 532208 : gel(z,2) = diviiexact(gel(y,1), p1);
2580 532208 : gel(z,1) = mulii(gel(y,2), x);
2581 532208 : normalize_frac(z);
2582 532208 : fix_frac_if_int_GC(z,tetpil);
2583 : }
2584 783374 : return z;
2585 :
2586 469 : case t_FFELT: return Z_FF_div(x,y);
2587 669857 : case t_COMPLEX: return divRc(x,y);
2588 1505 : case t_PADIC: return divTp(x, y);
2589 231 : case t_QUAD:
2590 231 : av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
2591 231 : return gerepile(av, tetpil, gdiv(p2,p1));
2592 : }
2593 : }
2594 217691488 : if (gequal0(y))
2595 : {
2596 49 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2597 28 : if (ty != t_MAT) pari_err_INV("gdiv",y);
2598 : }
2599 :
2600 217701254 : if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2601 : {
2602 181427379 : case t_REAL:
2603 181427379 : switch(ty)
2604 : {
2605 179060504 : case t_INT: return divri(x,y);
2606 519124 : case t_FRAC:
2607 519124 : av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2608 519124 : return gerepileuptoleaf(av, z);
2609 1847823 : case t_COMPLEX: return divRc(x, y);
2610 42 : case t_QUAD: return divfq(x, y, realprec(x));
2611 18 : default: pari_err_TYPE2("/",x,y);
2612 : }
2613 :
2614 301 : case t_INTMOD:
2615 : switch(ty)
2616 : {
2617 210 : case t_INT:
2618 210 : z = cgetg(3, t_INTMOD);
2619 210 : return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2620 35 : case t_FRAC: { GEN X = gel(x,1);
2621 35 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2622 35 : return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2623 : }
2624 7 : case t_FFELT:
2625 7 : if (!equalii(gel(x,1),FF_p_i(y)))
2626 0 : pari_err_OP("/",x,y);
2627 7 : return Z_FF_div(gel(x,2),y);
2628 :
2629 0 : case t_COMPLEX:
2630 0 : av = avma;
2631 0 : return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
2632 :
2633 14 : case t_QUAD:
2634 14 : av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
2635 14 : return gerepile(av,tetpil, gdiv(p2,p1));
2636 :
2637 7 : case t_PADIC: { GEN X = gel(x,1);
2638 7 : z = cgetg(3, t_INTMOD);
2639 7 : return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2640 : }
2641 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2642 : }
2643 :
2644 : case t_FRAC:
2645 : switch(ty)
2646 : {
2647 2027389 : case t_INT: z = cgetg(3, t_FRAC);
2648 2027389 : p1 = gcdii(y,gel(x,1));
2649 2027389 : if (equali1(p1))
2650 : {
2651 781973 : set_avma((pari_sp)z); tetpil = 0;
2652 781973 : gel(z,1) = icopy(gel(x,1));
2653 : }
2654 : else
2655 : {
2656 1245416 : y = diviiexact(y,p1); tetpil = avma;
2657 1245416 : gel(z,1) = diviiexact(gel(x,1), p1);
2658 : }
2659 2027389 : gel(z,2) = mulii(gel(x,2),y);
2660 2027389 : normalize_frac(z);
2661 2027389 : if (tetpil) fix_frac_if_int_GC(z,tetpil);
2662 2027389 : return z;
2663 :
2664 57742 : case t_REAL:
2665 57742 : av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
2666 57742 : return gerepile(av, tetpil, divir(gel(x,1), p1));
2667 :
2668 7 : case t_INTMOD: { GEN Y = gel(y,1);
2669 7 : z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2670 7 : return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2671 : }
2672 :
2673 28 : case t_FFELT: av=avma;
2674 28 : return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2675 :
2676 483 : case t_COMPLEX: return divRc(x, y);
2677 :
2678 2128 : case t_PADIC:
2679 2128 : if (!signe(gel(x,1))) return gen_0;
2680 2128 : return divTp(x, y);
2681 :
2682 42 : case t_QUAD:
2683 42 : av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
2684 42 : return gerepile(av,tetpil,gdiv(p2,p1));
2685 : }
2686 :
2687 : case t_FFELT:
2688 133 : switch (ty)
2689 : {
2690 49 : case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2691 28 : case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2692 7 : case t_INTMOD:
2693 7 : if (!equalii(gel(y,1),FF_p_i(x)))
2694 0 : pari_err_OP("/",x,y);
2695 7 : return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2696 49 : default:
2697 49 : pari_err_TYPE2("/",x,y);
2698 : }
2699 0 : break;
2700 :
2701 12477848 : case t_COMPLEX:
2702 : switch(ty)
2703 : {
2704 12477825 : case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2705 0 : case t_INTMOD: return mulRc_direct(ginv(y), x);
2706 0 : case t_PADIC:
2707 0 : return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
2708 0 : case t_QUAD:
2709 0 : lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2710 0 : return divfq(x, y, lx);
2711 : }
2712 :
2713 : case t_PADIC:
2714 : switch(ty)
2715 : {
2716 1327641 : case t_INT: case t_FRAC: { GEN p = gel(x,2);
2717 1327543 : return signe(gel(x,4))? divpT(x, y)
2718 2655184 : : zeropadic(p, valp(x) - Q_pval(y,p));
2719 : }
2720 7 : case t_INTMOD: { GEN Y = gel(y,1);
2721 7 : z = cgetg(3, t_INTMOD);
2722 7 : return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2723 : }
2724 14 : case t_COMPLEX: case t_QUAD:
2725 14 : av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
2726 14 : return gerepile(av,tetpil,gdiv(p1,p2));
2727 :
2728 28 : case t_REAL: pari_err_TYPE2("/",x,y);
2729 : }
2730 :
2731 : case t_QUAD:
2732 : switch (ty)
2733 : {
2734 1190 : case t_INT: case t_INTMOD: case t_FRAC:
2735 1190 : z = cgetg(4,t_QUAD);
2736 1190 : gel(z,1) = ZX_copy(gel(x,1));
2737 1190 : gel(z,2) = gdiv(gel(x,2), y);
2738 1190 : gel(z,3) = gdiv(gel(x,3), y); return z;
2739 63 : case t_REAL: return divqf(x, y, realprec(y));
2740 14 : case t_PADIC: return divTp(x, y);
2741 0 : case t_COMPLEX:
2742 0 : ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2743 0 : return divqf(x, y, ly);
2744 : }
2745 : }
2746 20378615 : switch(ty) {
2747 377327 : case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2748 377327 : return gmul(x, ginv(y)); /* missing gerepile, for speed */
2749 42 : case t_MAT:
2750 42 : av = avma; p1 = RgM_inv(y);
2751 35 : if (!p1) pari_err_INV("gdiv",y);
2752 28 : return gerepileupto(av, gmul(x, p1));
2753 0 : case t_VEC: case t_COL:
2754 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2755 0 : pari_err_TYPE2("/",x,y);
2756 : }
2757 20001246 : switch(tx) {
2758 3197288 : case t_VEC: case t_COL: case t_MAT:
2759 3197288 : z = cgetg_copy(x, &lx);
2760 11186514 : for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
2761 3197264 : return z;
2762 0 : case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2763 0 : pari_err_TYPE2("/",x,y);
2764 : }
2765 :
2766 16803958 : vy = gvar(y);
2767 16804724 : if (tx == t_POLMOD) { GEN X = gel(x,1);
2768 207982 : vx = varn(X);
2769 207982 : if (vx != vy) {
2770 207408 : if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2771 207079 : retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2772 : }
2773 : /* y is POL, SER or RFRAC */
2774 574 : av = avma;
2775 574 : switch(ty)
2776 : {
2777 0 : case t_RFRAC: y = gmod(ginv(y), X); break;
2778 574 : default: y = ginvmod(gmod(y,X), X);
2779 : }
2780 567 : return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
2781 : }
2782 : /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2783 : * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2784 : * variable NO_VARIABLE, then the operation is incorrect */
2785 16596742 : vx = gvar(x);
2786 16596748 : if (vx != vy) { /* includes cases where one is scalar */
2787 14701598 : if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2788 9267200 : else return div_scal_T(x, y, ty);
2789 : }
2790 1895150 : switch(tx)
2791 : {
2792 1046164 : case t_POL:
2793 : switch(ty)
2794 : {
2795 28 : case t_SER:
2796 : {
2797 28 : GEN y0 = y;
2798 : long v;
2799 28 : av = avma; v = RgX_valrem(x, &x);
2800 28 : if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
2801 14 : v -= valser(y); ly = lg(y); /* > 2 */
2802 14 : y = ser2pol_i_normalize(y, ly, &i);
2803 14 : if (i)
2804 : {
2805 7 : pari_warn(warner,"normalizing a series with 0 leading term");
2806 7 : ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
2807 : }
2808 7 : z = init_ser(ly, vx, v);
2809 7 : return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, ly-2)));
2810 : }
2811 :
2812 1046136 : case t_RFRAC:
2813 : {
2814 1046136 : GEN y1 = gel(y,1), y2 = gel(y,2);
2815 1046136 : if (typ(y1) == t_POL && varn(y1) == vx)
2816 35 : return mul_rfrac_scal(y2, y1, x);
2817 1046101 : av = avma;
2818 1046101 : return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2819 : }
2820 : }
2821 0 : break;
2822 :
2823 578143 : case t_SER:
2824 : switch(ty)
2825 : {
2826 578129 : case t_POL:
2827 : {
2828 578129 : long v = valser(x);
2829 578129 : lx = lg(x);
2830 578129 : if (lx == 2) return zeroser(vx, v - RgX_val(y));
2831 578024 : av = avma;
2832 578024 : x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
2833 578024 : z = init_ser(lx, vx, v);
2834 578024 : if (!signe(x)) setsigne(z,0);
2835 578024 : return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
2836 : }
2837 14 : case t_RFRAC:
2838 14 : av = avma;
2839 14 : return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2840 : }
2841 0 : break;
2842 :
2843 270821 : case t_RFRAC:
2844 : switch(ty)
2845 : {
2846 270821 : case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2847 0 : case t_SER:
2848 0 : av = avma;
2849 0 : return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
2850 : }
2851 0 : break;
2852 : }
2853 22 : pari_err_TYPE2("/",x,y);
2854 : return NULL; /* LCOV_EXCL_LINE */
2855 : }
2856 :
2857 : /********************************************************************/
2858 : /** **/
2859 : /** SIMPLE MULTIPLICATION **/
2860 : /** **/
2861 : /********************************************************************/
2862 : GEN
2863 248608517 : gmulsg(long s, GEN y)
2864 : {
2865 : long ly, i;
2866 : pari_sp av;
2867 : GEN z;
2868 :
2869 248608517 : switch(typ(y))
2870 : {
2871 157074670 : case t_INT: return mulsi(s,y);
2872 70287151 : case t_REAL: return s? mulsr(s,y): gen_0; /* gmul semantic */
2873 174336 : case t_INTMOD: { GEN p = gel(y,1);
2874 174336 : z = cgetg(3,t_INTMOD);
2875 174337 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
2876 174335 : gel(z,1) = icopy(p); return z;
2877 : }
2878 548265 : case t_FFELT: return FF_Z_mul(y,stoi(s));
2879 5713766 : case t_FRAC:
2880 5713766 : if (!s) return gen_0;
2881 5708810 : z = cgetg(3,t_FRAC);
2882 5708810 : i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
2883 5708810 : if (i == 1)
2884 : {
2885 2184727 : gel(z,2) = icopy(gel(y,2));
2886 2184727 : gel(z,1) = mulis(gel(y,1), s);
2887 : }
2888 : else
2889 : {
2890 3524083 : gel(z,2) = diviuexact(gel(y,2), (ulong)i);
2891 3524083 : gel(z,1) = mulis(gel(y,1), s/i);
2892 3524083 : fix_frac_if_int(z);
2893 : }
2894 5708810 : return z;
2895 :
2896 13317123 : case t_COMPLEX:
2897 13317123 : if (!s) return gen_0;
2898 13242775 : z = cgetg(3, t_COMPLEX);
2899 13242674 : gel(z,1) = gmulsg(s,gel(y,1));
2900 13241127 : gel(z,2) = gmulsg(s,gel(y,2)); return z;
2901 :
2902 1435 : case t_PADIC:
2903 1435 : if (!s) return gen_0;
2904 1435 : av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
2905 :
2906 7 : case t_QUAD: z = cgetg(4, t_QUAD);
2907 7 : gel(z,1) = ZX_copy(gel(y,1));
2908 7 : gel(z,2) = gmulsg(s,gel(y,2));
2909 7 : gel(z,3) = gmulsg(s,gel(y,3)); return z;
2910 :
2911 75173 : case t_POLMOD:
2912 75173 : retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
2913 :
2914 594196 : case t_POL:
2915 594196 : if (!signe(y)) return RgX_copy(y);
2916 579629 : if (!s) return scalarpol(Rg_get_0(y), varn(y));
2917 577984 : z = cgetg_copy(y, &ly); z[1]=y[1];
2918 2387520 : for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2919 577984 : return normalizepol_lg(z, ly);
2920 :
2921 182 : case t_SER:
2922 182 : if (ser_isexactzero(y)) return gcopy(y);
2923 182 : if (!s) return Rg_get_0(y);
2924 182 : z = cgetg_copy(y, &ly); z[1]=y[1];
2925 3864 : for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2926 182 : return normalizeser(z);
2927 :
2928 14 : case t_RFRAC:
2929 14 : if (!s) return zeropol(varn(gel(y,2)));
2930 14 : if (s == 1) return gcopy(y);
2931 14 : if (s == -1) return gneg(y);
2932 14 : return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
2933 :
2934 832396 : case t_VEC: case t_COL: case t_MAT:
2935 832396 : z = cgetg_copy(y, &ly);
2936 2273907 : for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2937 832391 : return z;
2938 : }
2939 0 : pari_err_TYPE("gmulsg",y);
2940 : return NULL; /* LCOV_EXCL_LINE */
2941 : }
2942 :
2943 : GEN
2944 133754682 : gmulug(ulong s, GEN y)
2945 : {
2946 : long ly, i;
2947 : pari_sp av;
2948 : GEN z;
2949 :
2950 133754682 : switch(typ(y))
2951 : {
2952 131703394 : case t_INT: return mului(s,y);
2953 1054628 : case t_REAL: return s? mulur(s,y): gen_0; /* gmul semantic */
2954 364 : case t_INTMOD: { GEN p = gel(y,1);
2955 364 : z = cgetg(3,t_INTMOD);
2956 364 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(mului(s,gel(y,2)), p));
2957 364 : gel(z,1) = icopy(p); return z;
2958 : }
2959 413 : case t_FFELT: return FF_Z_mul(y,utoi(s));
2960 937079 : case t_FRAC:
2961 937079 : if (!s) return gen_0;
2962 937065 : z = cgetg(3,t_FRAC);
2963 937065 : i = ugcd(s, umodiu(gel(y,2), s));
2964 937065 : if (i == 1)
2965 : {
2966 793453 : gel(z,2) = icopy(gel(y,2));
2967 793453 : gel(z,1) = muliu(gel(y,1), s);
2968 : }
2969 : else
2970 : {
2971 143612 : gel(z,2) = diviuexact(gel(y,2), i);
2972 143612 : gel(z,1) = muliu(gel(y,1), s/i);
2973 143612 : fix_frac_if_int(z);
2974 : }
2975 937065 : return z;
2976 :
2977 24799 : case t_COMPLEX:
2978 24799 : if (!s) return gen_0;
2979 24799 : z = cgetg(3, t_COMPLEX);
2980 24799 : gel(z,1) = gmulug(s,gel(y,1));
2981 24799 : gel(z,2) = gmulug(s,gel(y,2)); return z;
2982 :
2983 3898 : case t_PADIC:
2984 3898 : if (!s) return gen_0;
2985 3898 : av = avma; return gerepileupto(av, mulpp(cvtop2(utoi(s),y), y));
2986 :
2987 0 : case t_QUAD: z = cgetg(4, t_QUAD);
2988 0 : gel(z,1) = ZX_copy(gel(y,1));
2989 0 : gel(z,2) = gmulug(s,gel(y,2));
2990 0 : gel(z,3) = gmulug(s,gel(y,3)); return z;
2991 :
2992 3199 : case t_POLMOD:
2993 3199 : retmkpolmod(gmulug(s,gel(y,2)), RgX_copy(gel(y,1)));
2994 :
2995 12642 : case t_POL:
2996 12642 : if (!signe(y)) return RgX_copy(y);
2997 11865 : if (!s) return scalarpol(Rg_get_0(y), varn(y));
2998 11865 : z = cgetg_copy(y, &ly); z[1]=y[1];
2999 34923 : for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
3000 11865 : return normalizepol_lg(z, ly);
3001 :
3002 0 : case t_SER:
3003 0 : if (ser_isexactzero(y)) return gcopy(y);
3004 0 : if (!s) return Rg_get_0(y);
3005 0 : z = cgetg_copy(y, &ly); z[1]=y[1];
3006 0 : for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
3007 0 : return normalizeser(z);
3008 :
3009 0 : case t_RFRAC:
3010 0 : if (!s) return zeropol(varn(gel(y,2)));
3011 0 : if (s == 1) return gcopy(y);
3012 0 : return mul_rfrac_scal(gel(y,1), gel(y,2), utoi(s));
3013 :
3014 14266 : case t_VEC: case t_COL: case t_MAT:
3015 14266 : z = cgetg_copy(y, &ly);
3016 74284 : for (i=1; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
3017 14266 : return z;
3018 : }
3019 0 : pari_err_TYPE("gmulsg",y);
3020 : return NULL; /* LCOV_EXCL_LINE */
3021 : }
3022 :
3023 : /********************************************************************/
3024 : /** **/
3025 : /** SIMPLE DIVISION **/
3026 : /** **/
3027 : /********************************************************************/
3028 :
3029 : GEN
3030 8644011 : gdivgs(GEN x, long s)
3031 : {
3032 8644011 : long tx = typ(x), lx, i;
3033 : pari_sp av;
3034 : GEN z;
3035 :
3036 8644011 : if (!s)
3037 : {
3038 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3039 0 : pari_err_INV("gdivgs",gen_0);
3040 : }
3041 8644019 : switch(tx)
3042 : {
3043 1434225 : case t_INT: return Qdivis(x, s);
3044 5888575 : case t_REAL: return divrs(x,s);
3045 :
3046 357 : case t_INTMOD:
3047 357 : z = cgetg(3, t_INTMOD);
3048 357 : return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
3049 :
3050 735 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
3051 :
3052 498734 : case t_FRAC: z = cgetg(3, t_FRAC);
3053 498734 : i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
3054 498734 : if (i == 1)
3055 : {
3056 374017 : gel(z,2) = mulsi(s, gel(x,2));
3057 374017 : gel(z,1) = icopy(gel(x,1));
3058 : }
3059 : else
3060 : {
3061 124717 : gel(z,2) = mulsi(s/i, gel(x,2));
3062 124717 : gel(z,1) = divis(gel(x,1), i);
3063 : }
3064 498734 : normalize_frac(z);
3065 498734 : fix_frac_if_int(z); return z;
3066 :
3067 760171 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3068 760172 : gel(z,1) = gdivgs(gel(x,1),s);
3069 760168 : gel(z,2) = gdivgs(gel(x,2),s); return z;
3070 :
3071 133 : case t_PADIC: /* divpT */
3072 : {
3073 133 : GEN p = gel(x,2);
3074 133 : if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
3075 133 : av = avma;
3076 133 : return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
3077 : }
3078 :
3079 28 : case t_QUAD: z = cgetg(4, t_QUAD);
3080 28 : gel(z,1) = ZX_copy(gel(x,1));
3081 28 : gel(z,2) = gdivgs(gel(x,2),s);
3082 28 : gel(z,3) = gdivgs(gel(x,3),s); return z;
3083 :
3084 29666 : case t_POLMOD:
3085 29666 : retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
3086 :
3087 91 : case t_RFRAC:
3088 91 : if (s == 1) return gcopy(x);
3089 84 : else if (s == -1) return gneg(x);
3090 84 : return div_rfrac_scal(x, stoi(s));
3091 :
3092 29463 : case t_POL: case t_SER:
3093 29463 : z = cgetg_copy(x, &lx); z[1] = x[1];
3094 117166 : for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
3095 29463 : return z;
3096 1841 : case t_VEC: case t_COL: case t_MAT:
3097 1841 : z = cgetg_copy(x, &lx);
3098 41335 : for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
3099 1841 : return z;
3100 :
3101 : }
3102 0 : pari_err_TYPE2("/",x, stoi(s));
3103 : return NULL; /* LCOV_EXCL_LINE */
3104 : }
3105 :
3106 : GEN
3107 43481882 : gdivgu(GEN x, ulong s)
3108 : {
3109 43481882 : long tx = typ(x), lx, i;
3110 : pari_sp av;
3111 : GEN z;
3112 :
3113 43481882 : if (!s)
3114 : {
3115 0 : if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
3116 0 : pari_err_INV("gdivgu",gen_0);
3117 : }
3118 43481905 : switch(tx)
3119 : {
3120 17274519 : case t_INT: return Qdiviu(x, s);
3121 12895703 : case t_REAL: return divru(x,s);
3122 :
3123 210315 : case t_INTMOD:
3124 210315 : z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
3125 210315 : return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
3126 :
3127 308 : case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
3128 :
3129 637355 : case t_FRAC: z = cgetg(3, t_FRAC);
3130 637355 : i = ugcd(s, umodiu(gel(x,1), s));
3131 637355 : if (i == 1)
3132 : {
3133 457990 : gel(z,2) = mului(s, gel(x,2));
3134 457990 : gel(z,1) = icopy(gel(x,1));
3135 : }
3136 : else
3137 : {
3138 179365 : gel(z,2) = mului(s/i, gel(x,2));
3139 179365 : gel(z,1) = divis(gel(x,1), i);
3140 : }
3141 637355 : normalize_frac(z);
3142 637355 : fix_frac_if_int(z); return z;
3143 :
3144 11557599 : case t_COMPLEX: z = cgetg(3, t_COMPLEX);
3145 11557599 : gel(z,1) = gdivgu(gel(x,1),s);
3146 11557598 : gel(z,2) = gdivgu(gel(x,2),s); return z;
3147 :
3148 851940 : case t_PADIC: /* divpT */
3149 : {
3150 851940 : GEN p = gel(x,2);
3151 851940 : if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
3152 717834 : av = avma;
3153 717834 : return gerepileupto(av, divpp(x, cvtop2(utoi(s),x)));
3154 : }
3155 :
3156 0 : case t_QUAD: z = cgetg(4, t_QUAD);
3157 0 : gel(z,1) = ZX_copy(gel(x,1));
3158 0 : gel(z,2) = gdivgu(gel(x,2),s);
3159 0 : gel(z,3) = gdivgu(gel(x,3),s); return z;
3160 :
3161 1456 : case t_POLMOD:
3162 1456 : retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
3163 :
3164 56 : case t_RFRAC:
3165 56 : if (s == 1) return gcopy(x);
3166 56 : return div_rfrac_scal(x, utoi(s));
3167 :
3168 52353 : case t_POL: case t_SER:
3169 52353 : z = cgetg_copy(x, &lx); z[1] = x[1];
3170 230244 : for (i=2; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
3171 52353 : return z;
3172 301 : case t_VEC: case t_COL: case t_MAT:
3173 301 : z = cgetg_copy(x, &lx);
3174 1092 : for (i=1; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
3175 301 : return z;
3176 : }
3177 0 : pari_err_TYPE2("/",x, utoi(s));
3178 : return NULL; /* LCOV_EXCL_LINE */
3179 : }
3180 :
3181 : /* x / (i*(i+1)) */
3182 : GEN
3183 106508918 : divrunextu(GEN x, ulong i)
3184 : {
3185 106508918 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3186 0 : return divri(x, muluu(i , i+1));
3187 : else
3188 106508918 : return divru(x, i*(i+1));
3189 : }
3190 : /* x / (i*(i+1)) */
3191 : GEN
3192 804801 : gdivgunextu(GEN x, ulong i)
3193 : {
3194 804801 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
3195 0 : return gdivgu(x, i*(i+1));
3196 : else
3197 804801 : return gdiv(x, muluu(i, i+1));
3198 : }
3199 :
3200 : /* True shift (exact multiplication by 2^n) */
3201 : GEN
3202 89877037 : gmul2n(GEN x, long n)
3203 : {
3204 : long lx, i, k, l;
3205 : GEN z, a, b;
3206 :
3207 89877037 : switch(typ(x))
3208 : {
3209 20413120 : case t_INT:
3210 20413120 : if (n>=0) return shifti(x,n);
3211 3389568 : if (!signe(x)) return gen_0;
3212 3294246 : l = vali(x); n = -n;
3213 3294304 : if (n<=l) return shifti(x,-n);
3214 375945 : z = cgetg(3,t_FRAC);
3215 375945 : gel(z,1) = shifti(x,-l);
3216 375945 : gel(z,2) = int2n(n-l); return z;
3217 :
3218 50174243 : case t_REAL:
3219 50174243 : return shiftr(x,n);
3220 :
3221 181070 : case t_INTMOD: b = gel(x,1); a = gel(x,2);
3222 181070 : z = cgetg(3,t_INTMOD);
3223 181069 : if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3224 82908 : gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
3225 82907 : gel(z,1) = icopy(b); return z;
3226 :
3227 217440 : case t_FFELT: return FF_mul2n(x,n);
3228 :
3229 842911 : case t_FRAC: a = gel(x,1); b = gel(x,2);
3230 842911 : l = vali(a);
3231 842911 : k = vali(b);
3232 842911 : if (n+l >= k)
3233 : {
3234 301733 : if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3235 239642 : l = n-k; k = -k;
3236 : }
3237 : else
3238 : {
3239 541178 : k = -(l+n); l = -l;
3240 : }
3241 780820 : z = cgetg(3,t_FRAC);
3242 780820 : gel(z,1) = shifti(a,l);
3243 780820 : gel(z,2) = shifti(b,k); return z;
3244 :
3245 11799609 : case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3246 11799605 : gel(z,1) = gmul2n(gel(x,1),n);
3247 11799623 : gel(z,2) = gmul2n(gel(x,2),n); return z;
3248 :
3249 105 : case t_QUAD: z = cgetg(4,t_QUAD);
3250 105 : gel(z,1) = ZX_copy(gel(x,1));
3251 105 : gel(z,2) = gmul2n(gel(x,2),n);
3252 105 : gel(z,3) = gmul2n(gel(x,3),n); return z;
3253 :
3254 44219 : case t_POLMOD:
3255 44219 : retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3256 :
3257 1581299 : case t_POL:
3258 1581299 : z = cgetg_copy(x, &lx); z[1] = x[1];
3259 8600661 : for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3260 1581201 : return normalizepol_lg(z, lx); /* needed if char = 2 */
3261 98899 : case t_SER:
3262 98899 : if (ser_isexactzero(x)) return gcopy(x);
3263 98871 : z = cgetg_copy(x, &lx); z[1] = x[1];
3264 619311 : for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3265 98871 : return normalizeser(z); /* needed if char = 2 */
3266 4521046 : case t_VEC: case t_COL: case t_MAT:
3267 4521046 : z = cgetg_copy(x, &lx);
3268 17718509 : for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3269 4520999 : return z;
3270 :
3271 21 : case t_RFRAC: /* int2n wrong if n < 0 */
3272 21 : return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3273 :
3274 4536 : case t_PADIC: /* int2n wrong if n < 0 */
3275 4536 : return gmul(gmul2n(gen_1,n),x);
3276 : }
3277 0 : pari_err_TYPE("gmul2n",x);
3278 : return NULL; /* LCOV_EXCL_LINE */
3279 : }
3280 :
3281 : /*******************************************************************/
3282 : /* */
3283 : /* INVERSE */
3284 : /* */
3285 : /*******************************************************************/
3286 : static GEN
3287 79929 : inv_polmod(GEN T, GEN x)
3288 : {
3289 79929 : GEN z = cgetg(3,t_POLMOD), a;
3290 79928 : gel(z,1) = RgX_copy(T);
3291 79930 : if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3292 62650 : a = ginv(x);
3293 : else
3294 : {
3295 17280 : if (lg(T) == 5) /* quadratic fields */
3296 13037 : a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3297 : else
3298 4243 : a = RgXQ_inv(x, T);
3299 : }
3300 79929 : gel(z,2) = a; return z;
3301 : }
3302 :
3303 : GEN
3304 23298208 : ginv(GEN x)
3305 : {
3306 : long s;
3307 : pari_sp av, tetpil;
3308 : GEN z, y, p1, p2;
3309 :
3310 23298208 : switch(typ(x))
3311 : {
3312 4110313 : case t_INT:
3313 4110313 : if (is_pm1(x)) return icopy(x);
3314 2582539 : s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3315 2582525 : z = cgetg(3,t_FRAC);
3316 2582560 : gel(z,1) = s<0? gen_m1: gen_1;
3317 2582560 : gel(z,2) = absi(x); return z;
3318 :
3319 5955523 : case t_REAL: return invr(x);
3320 :
3321 7035 : case t_INTMOD: z=cgetg(3,t_INTMOD);
3322 7035 : gel(z,1) = icopy(gel(x,1));
3323 7035 : gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3324 :
3325 425856 : case t_FRAC: {
3326 425856 : GEN a = gel(x,1), b = gel(x,2);
3327 425856 : s = signe(a);
3328 425856 : if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3329 243331 : z = cgetg(3,t_FRAC);
3330 243331 : gel(z,1) = icopy(b);
3331 243331 : gel(z,2) = icopy(a);
3332 243331 : normalize_frac(z); return z;
3333 : }
3334 6881470 : case t_COMPLEX:
3335 6881470 : av=avma;
3336 6881470 : p1=cxnorm(x);
3337 6880617 : p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
3338 6881301 : tetpil=avma;
3339 6881301 : return gerepile(av,tetpil,divcR(p2,p1));
3340 :
3341 273 : case t_QUAD:
3342 273 : av=avma; p1=quadnorm(x); p2=conj_i(x); tetpil=avma;
3343 273 : return gerepile(av,tetpil,gdiv(p2,p1));
3344 :
3345 346227 : case t_PADIC: z = cgetg(5,t_PADIC);
3346 346227 : if (!signe(gel(x,4))) pari_err_INV("ginv",x);
3347 346220 : z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
3348 346220 : gel(z,2) = icopy(gel(x,2));
3349 346220 : gel(z,3) = icopy(gel(x,3));
3350 346220 : gel(z,4) = Zp_inv(gel(x,4),gel(z,2),precp(x)); return z;
3351 :
3352 79929 : case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3353 14272 : case t_FFELT: return FF_inv(x);
3354 5321009 : case t_POL: return gred_rfrac_simple(gen_1,x);
3355 19418 : case t_SER: return ser_inv(x);
3356 2926 : case t_RFRAC:
3357 : {
3358 2926 : GEN n = gel(x,1), d = gel(x,2);
3359 2926 : pari_sp av = avma, ltop;
3360 2926 : if (gequal0(n)) pari_err_INV("ginv",x);
3361 :
3362 2926 : n = simplify_shallow(n);
3363 2926 : if (typ(n) != t_POL || varn(n) != varn(d))
3364 : {
3365 2926 : if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3366 679 : ltop = avma;
3367 679 : z = RgX_Rg_div(d,n);
3368 : } else {
3369 0 : ltop = avma;
3370 0 : z = cgetg(3,t_RFRAC);
3371 0 : gel(z,1) = RgX_copy(d);
3372 0 : gel(z,2) = RgX_copy(n);
3373 : }
3374 679 : stackdummy(av, ltop);
3375 679 : return z;
3376 : }
3377 :
3378 21 : case t_VEC: if (!is_ext_qfr(x)) break;
3379 : case t_QFB:
3380 99415 : return qfbpow(x, gen_m1);
3381 34643 : case t_MAT:
3382 34643 : y = RgM_inv(x);
3383 34636 : if (!y) pari_err_INV("ginv",x);
3384 34566 : return y;
3385 28 : case t_VECSMALL:
3386 : {
3387 28 : long i, lx = lg(x)-1;
3388 28 : y = zero_zv(lx);
3389 112 : for (i=1; i<=lx; i++)
3390 : {
3391 84 : long xi = x[i];
3392 84 : if (xi<1 || xi>lx || y[xi])
3393 0 : pari_err_TYPE("ginv [not a permutation]", x);
3394 84 : y[xi] = i;
3395 : }
3396 28 : return y;
3397 : }
3398 : }
3399 6 : pari_err_TYPE("inverse",x);
3400 : return NULL; /* LCOV_EXCL_LINE */
3401 : }
|
