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 /* assume z[1] was created last */
25 #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
26 z = gerepileupto((pari_sp)(z+3), gel(z,1));
27
28 /* assume z[1] was created last */
29 #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
30 z = gerepileupto((pari_sp)(z+3), gel(z,1));\
31 else\
32 gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
33
34 static void
warn_coercion(GEN x,GEN y,GEN z)35 warn_coercion(GEN x, GEN y, GEN z)
36 {
37 if (DEBUGLEVEL)
38 pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
39 }
40
41 static long
kro_quad(GEN x,GEN y)42 kro_quad(GEN x, GEN y)
43 { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
44
45 /* is -1 not a square in Zp, assume p prime */
46 INLINE int
Zp_nosquare_m1(GEN p)47 Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
48
49 static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
50 static GEN mulpp(GEN x, GEN y);
51 static GEN divpp(GEN x, GEN y);
52 /* Argument codes for inline routines
53 * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
54 * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
55 * T: some type (to be converted to PADIC)
56 */
57 static GEN
addRc(GEN x,GEN y)58 addRc(GEN x, GEN y) {
59 GEN z = cgetg(3,t_COMPLEX);
60 gel(z,1) = gadd(x,gel(y,1));
61 gel(z,2) = gcopy(gel(y,2)); return z;
62 }
63 static GEN
mulRc(GEN x,GEN y)64 mulRc(GEN x, GEN y) {
65 GEN z = cgetg(3,t_COMPLEX);
66 gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
67 gel(z,2) = gmul(x,gel(y,2)); return z;
68 }
69 /* for INTMODs: can't simplify when Re(x) = gen_0 */
70 static GEN
mulRc_direct(GEN x,GEN y)71 mulRc_direct(GEN x, GEN y) {
72 GEN z = cgetg(3,t_COMPLEX);
73 gel(z,1) = gmul(x,gel(y,1));
74 gel(z,2) = gmul(x,gel(y,2)); return z;
75 }
76 static GEN
divRc(GEN x,GEN y)77 divRc(GEN x, GEN y) {
78 GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
79 GEN z = cgetg(3,t_COMPLEX);
80 gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
81 gel(z,2) = gmul(mt, gel(y,2));
82 return z;
83 }
84 static GEN
divcR(GEN x,GEN y)85 divcR(GEN x, GEN y) {
86 GEN z = cgetg(3,t_COMPLEX);
87 gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
88 gel(z,2) = gdiv(gel(x,2), y); return z;
89 }
90 static GEN
addRq(GEN x,GEN y)91 addRq(GEN x, GEN y) {
92 GEN z = cgetg(4,t_QUAD);
93 gel(z,1) = ZX_copy(gel(y,1));
94 gel(z,2) = gadd(x, gel(y,2));
95 gel(z,3) = gcopy(gel(y,3)); return z;
96 }
97 static GEN
mulRq(GEN x,GEN y)98 mulRq(GEN x, GEN y) {
99 GEN z = cgetg(4,t_QUAD);
100 gel(z,1) = ZX_copy(gel(y,1));
101 gel(z,2) = gmul(x,gel(y,2));
102 gel(z,3) = gmul(x,gel(y,3)); return z;
103 }
104 static GEN
addqf(GEN x,GEN y,long prec)105 addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
106 long i = gexpo(x) - gexpo(y);
107 if (i > 0) prec += nbits2extraprec( i );
108 return gerepileupto(av, gadd(y, quadtofp(x, prec)));
109 }
110 static GEN
mulrfrac(GEN x,GEN y)111 mulrfrac(GEN x, GEN y)
112 {
113 pari_sp av;
114 GEN z, a = gel(y,1), b = gel(y,2);
115 if (is_pm1(a)) /* frequent special case */
116 {
117 z = divri(x, b); if (signe(a) < 0) togglesign(z);
118 return z;
119 }
120 av = avma;
121 return gerepileuptoleaf(av, divri(mulri(x,a), b));
122 }
123 static GEN
mulqf(GEN x,GEN y,long prec)124 mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
125 return gerepileupto(av, gmul(y, quadtofp(x, prec)));
126 }
127 static GEN
divqf(GEN x,GEN y,long prec)128 divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
129 return gerepileupto(av, gdiv(quadtofp(x,prec), y));
130 }
131 static GEN
divfq(GEN x,GEN y,long prec)132 divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
133 return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
134 }
135 /* y PADIC, x + y by converting x to padic */
136 static GEN
addTp(GEN x,GEN y)137 addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
138
139 if (!valp(y)) z = cvtop2(x,y);
140 else {
141 long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
142 z = cvtop(x, gel(y,2), l);
143 }
144 return gerepileupto(av, addsub_pp(z, y, addii));
145 }
146 /* y PADIC, x * y by converting x to padic */
147 static GEN
mulTp(GEN x,GEN y)148 mulTp(GEN x, GEN y) { pari_sp av = avma;
149 return gerepileupto(av, mulpp(cvtop2(x,y), y));
150 }
151 /* y PADIC, non zero x / y by converting x to padic */
152 static GEN
divTp(GEN x,GEN y)153 divTp(GEN x, GEN y) { pari_sp av = avma;
154 return gerepileupto(av, divpp(cvtop2(x,y), y));
155 }
156 /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
157 * converted to O(p^e) and division by 0 */
158 static GEN
divpT(GEN x,GEN y)159 divpT(GEN x, GEN y) { pari_sp av = avma;
160 return gerepileupto(av, divpp(x, cvtop2(y,x)));
161 }
162
163 /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
164 * clean memory from z on */
165 static GEN
add_intmod_same(GEN z,GEN X,GEN x,GEN y)166 add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
167 if (lgefint(X) == 3) {
168 ulong u = Fl_add(itou(x),itou(y), X[2]);
169 set_avma((pari_sp)z); gel(z,2) = utoi(u);
170 }
171 else {
172 GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
173 gel(z,2) = gerepileuptoint((pari_sp)z, u);
174 }
175 gel(z,1) = icopy(X); return z;
176 }
177 static GEN
sub_intmod_same(GEN z,GEN X,GEN x,GEN y)178 sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
179 if (lgefint(X) == 3) {
180 ulong u = Fl_sub(itou(x),itou(y), X[2]);
181 set_avma((pari_sp)z); gel(z,2) = utoi(u);
182 }
183 else {
184 GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
185 gel(z,2) = gerepileuptoint((pari_sp)z, u);
186 }
187 gel(z,1) = icopy(X); return z;
188 }
189 /* cf add_intmod_same */
190 static GEN
mul_intmod_same(GEN z,GEN X,GEN x,GEN y)191 mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
192 if (lgefint(X) == 3) {
193 ulong u = Fl_mul(itou(x),itou(y), X[2]);
194 set_avma((pari_sp)z); gel(z,2) = utoi(u);
195 }
196 else
197 gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
198 gel(z,1) = icopy(X); return z;
199 }
200 /* cf add_intmod_same */
201 static GEN
div_intmod_same(GEN z,GEN X,GEN x,GEN y)202 div_intmod_same(GEN z, GEN X, GEN x, GEN y)
203 {
204 if (lgefint(X) == 3) {
205 ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
206 set_avma((pari_sp)z); gel(z,2) = utoi(u);
207 }
208 else
209 gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
210 gel(z,1) = icopy(X); return z;
211 }
212
213 /*******************************************************************/
214 /* */
215 /* REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC) */
216 /* */
217 /* (static routines are not memory clean, but OK for gerepileupto) */
218 /*******************************************************************/
219 /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
220 * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
221 static GEN
rfrac_denom_mul_scal(GEN d,GEN y)222 rfrac_denom_mul_scal(GEN d, GEN y)
223 {
224 GEN D = RgX_Rg_mul(d, y);
225 if (lg(D) != lg(d))
226 { /* try to generate a meaningful diagnostic */
227 D = gdiv(leading_coeff(d), y); /* should fail */
228 pari_err_INV("gred_rfrac", y); /* better than nothing */
229 }
230 return D;
231 }
232
233 /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
234 GEN
gred_rfrac_simple(GEN n,GEN d)235 gred_rfrac_simple(GEN n, GEN d)
236 {
237 GEN c, cn, cd, z;
238 long dd = degpol(d);
239
240 if (dd <= 0)
241 {
242 if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
243 n = gdiv(n, gel(d,2));
244 if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
245 return n;
246 }
247
248 cd = content(d);
249 while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
250 cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
251 if (!gequal1(cd)) {
252 d = RgX_Rg_div(d,cd);
253 if (!gequal1(cn))
254 {
255 if (gequal0(cn)) {
256 if (isexactzero(cn)) return scalarpol(cn, varn(d));
257 n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
258 c = gen_1;
259 } else {
260 n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
261 c = gdiv(cn,cd);
262 }
263 }
264 else
265 c = ginv(cd);
266 } else {
267 if (!gequal1(cn))
268 {
269 if (gequal0(cn)) {
270 if (isexactzero(cn)) return scalarpol(cn, varn(d));
271 c = gen_1;
272 } else {
273 n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
274 c = cn;
275 }
276 } else {
277 GEN y = cgetg(3,t_RFRAC);
278 gel(y,1) = gcopy(n);
279 gel(y,2) = RgX_copy(d); return y;
280 }
281 }
282
283 if (typ(c) == t_POL)
284 {
285 z = c;
286 do { z = content(z); } while (typ(z) == t_POL);
287 cd = denom_i(z);
288 cn = gmul(c, cd);
289 }
290 else
291 {
292 cn = numer_i(c);
293 cd = denom_i(c);
294 }
295 z = cgetg(3,t_RFRAC);
296 gel(z,1) = gmul(n, cn);
297 gel(z,2) = rfrac_denom_mul_scal(d, cd);
298 return z;
299 }
300
301 /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
302 static GEN
fix_rfrac(GEN x,long d)303 fix_rfrac(GEN x, long d)
304 {
305 GEN z, N, D;
306 if (!d || typ(x) == t_POL) return x;
307 z = cgetg(3, t_RFRAC);
308 N = gel(x,1);
309 D = gel(x,2);
310 if (d > 0) {
311 gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
312 : monomialcopy(N,d,varn(D));
313 gel(z, 2) = RgX_copy(D);
314 } else {
315 gel(z, 1) = gcopy(N);
316 gel(z, 2) = RgX_shift(D, -d);
317 }
318 return z;
319 }
320
321 /* assume d != 0 */
322 static GEN
gred_rfrac2(GEN n,GEN d)323 gred_rfrac2(GEN n, GEN d)
324 {
325 GEN y, z;
326 long v, vd, vn;
327
328 n = simplify_shallow(n);
329 if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
330 d = simplify_shallow(d);
331 if (typ(d) != t_POL) return gdiv(n,d);
332 vd = varn(d);
333 if (typ(n) != t_POL)
334 {
335 if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
336 if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
337 pari_err_BUG("gred_rfrac2 [incompatible variables]");
338 }
339 vn = varn(n);
340 if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
341 if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
342
343 /* now n and d are t_POLs in the same variable */
344 v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
345 if (!degpol(d))
346 {
347 n = RgX_Rg_div(n,gel(d,2));
348 return v? RgX_mulXn(n,v): n;
349 }
350
351 /* X does not divide gcd(n,d), deg(d) > 0 */
352 if (!isinexact(n) && !isinexact(d))
353 {
354 y = RgX_divrem(n, d, &z);
355 if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
356 z = RgX_gcd(d, z);
357 if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
358 }
359 return fix_rfrac(gred_rfrac_simple(n,d), v);
360 }
361
362 /* x,y t_INT, return x/y in reduced form */
363 GEN
Qdivii(GEN x,GEN y)364 Qdivii(GEN x, GEN y)
365 {
366 pari_sp av = avma;
367 GEN r, q;
368
369 if (lgefint(y) == 3)
370 {
371 q = Qdiviu(x, y[2]);
372 if (signe(y) > 0) return q;
373 if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
374 return q;
375 }
376 if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
377 if (equali1(x))
378 {
379 if (!signe(y)) pari_err_INV("gdiv",y);
380 retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
381 }
382 q = dvmdii(x,y,&r);
383 if (r == gen_0) return q; /* gen_0 intended */
384 r = gcdii(y, r);
385 if (lgefint(r) == 3)
386 {
387 ulong t = r[2];
388 set_avma(av);
389 if (t == 1) q = mkfraccopy(x,y);
390 else
391 {
392 q = cgetg(3,t_FRAC);
393 gel(q,1) = diviuexact(x,t);
394 gel(q,2) = diviuexact(y,t);
395 }
396 }
397 else
398 { /* rare: r and q left on stack for efficiency */
399 q = cgetg(3,t_FRAC);
400 gel(q,1) = diviiexact(x,r);
401 gel(q,2) = diviiexact(y,r);
402 }
403 normalize_frac(q); return q;
404 }
405
406 /* x t_INT, return x/y in reduced form */
407 GEN
Qdiviu(GEN x,ulong y)408 Qdiviu(GEN x, ulong y)
409 {
410 pari_sp av = avma;
411 ulong r, t;
412 GEN q;
413
414 if (y == 1) return icopy(x);
415 if (!y) pari_err_INV("gdiv",gen_0);
416 if (equali1(x)) retmkfrac(gen_1, utoipos(y));
417 q = absdiviu_rem(x,y,&r);
418 if (!r)
419 {
420 if (signe(x) < 0) togglesign(q);
421 return q;
422 }
423 t = ugcd(y, r); set_avma(av);
424 if (t == 1) retmkfrac(icopy(x), utoipos(y));
425 retmkfrac(diviuexact(x,t), utoipos(y / t));
426 }
427
428 /* x t_INT, return x/y in reduced form */
429 GEN
Qdivis(GEN x,long y)430 Qdivis(GEN x, long y)
431 {
432 pari_sp av = avma;
433 ulong r, t;
434 long s;
435 GEN q;
436
437 if (y > 0) return Qdiviu(x, y);
438 if (!y) pari_err_INV("gdiv",gen_0);
439 s = signe(x);
440 if (!s) return gen_0;
441 if (y < 0) { y = -y; s = -s; }
442 if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
443 if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
444 q = absdiviu_rem(x,y,&r);
445 if (!r)
446 {
447 if (s < 0) togglesign(q);
448 return q;
449 }
450 t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
451 if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
452 gel(q,1) = x; setsigne(x, s);
453 gel(q,2) = utoipos(y); return q;
454 }
455
456 /*******************************************************************/
457 /* */
458 /* CONJUGATION */
459 /* */
460 /*******************************************************************/
461 /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
462 static GEN
quad_polmod_conj(GEN x,GEN y)463 quad_polmod_conj(GEN x, GEN y)
464 {
465 GEN z, u, v, a, b;
466 if (typ(x) != t_POL) return gcopy(x);
467 if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
468 a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
469 b = gel(y,3); v = gel(x,2);
470 z = cgetg(4, t_POL); z[1] = x[1];
471 gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
472 gel(z,3) = gneg(u); return z;
473 }
474 static GEN
quad_polmod_norm(GEN x,GEN y)475 quad_polmod_norm(GEN x, GEN y)
476 {
477 GEN z, u, v, a, b, c;
478 if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
479 a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
480 b = gel(y,3); v = gel(x,2);
481 c = gel(y,2);
482 z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
483 if (!gequal1(a)) z = gdiv(z, a);
484 return gadd(z, gsqr(v));
485 }
486
487 GEN
conj_i(GEN x)488 conj_i(GEN x)
489 {
490 long lx, i;
491 GEN y;
492
493 switch(typ(x))
494 {
495 case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
496 return x;
497
498 case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
499
500 case t_QUAD:
501 y = cgetg(4,t_QUAD);
502 gel(y,1) = gel(x,1);
503 gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
504 : gadd(gel(x,2), gel(x,3));
505 gel(y,3) = gneg(gel(x,3)); return y;
506
507 case t_POL: case t_SER:
508 y = cgetg_copy(x, &lx); y[1] = x[1];
509 for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
510 return y;
511
512 case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
513 y = cgetg_copy(x, &lx);
514 for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
515 return y;
516
517 case t_POLMOD:
518 {
519 GEN X = gel(x,1);
520 long d = degpol(X);
521 if (d < 2) return x;
522 if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
523 }
524 }
525 pari_err_TYPE("gconj",x);
526 return NULL; /* LCOV_EXCL_LINE */
527 }
528 GEN
gconj(GEN x)529 gconj(GEN x)
530 { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
531
532 GEN
conjvec(GEN x,long prec)533 conjvec(GEN x,long prec)
534 {
535 long lx, s, i;
536 GEN z;
537
538 switch(typ(x))
539 {
540 case t_INT: case t_INTMOD: case t_FRAC:
541 return mkcolcopy(x);
542
543 case t_COMPLEX: case t_QUAD:
544 z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
545
546 case t_FFELT:
547 return FF_conjvec(x);
548
549 case t_VEC: case t_COL:
550 lx = lg(x); z = cgetg(lx,t_MAT);
551 if (lx == 1) return z;
552 gel(z,1) = conjvec(gel(x,1),prec);
553 s = lgcols(z);
554 for (i=2; i<lx; i++)
555 {
556 gel(z,i) = conjvec(gel(x,i),prec);
557 if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
558 }
559 break;
560
561 case t_POLMOD: {
562 GEN T = gel(x,1), r;
563 pari_sp av;
564
565 lx = lg(T);
566 if (lx <= 3) return cgetg(1,t_COL);
567 x = gel(x,2);
568 for (i=2; i<lx; i++)
569 {
570 GEN c = gel(T,i);
571 switch(typ(c)) {
572 case t_INTMOD: {
573 GEN p = gel(c,1);
574 pari_sp av;
575 if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
576 av = avma;
577 T = RgX_to_FpX(T,p);
578 x = RgX_to_FpX(x, p);
579 if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
580 z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
581 return gerepileupto(av, z);
582 }
583 case t_INT:
584 case t_FRAC: break;
585 default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
586 }
587 }
588 if (typ(x) != t_POL)
589 {
590 if (!is_rational_t(typ(x)))
591 pari_err_TYPE("conjvec [not a rational t_POL]",x);
592 retconst_col(lx-3, gcopy(x));
593 }
594 RgX_check_QX(x,"conjvec");
595 av = avma;
596 if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
597 r = cleanroots(T,prec);
598 z = cgetg(lx-2,t_COL);
599 for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
600 return gerepileupto(av, z);
601 }
602
603 default:
604 pari_err_TYPE("conjvec",x);
605 return NULL; /* LCOV_EXCL_LINE */
606 }
607 return z;
608 }
609
610 /********************************************************************/
611 /** **/
612 /** ADDITION **/
613 /** **/
614 /********************************************************************/
615 /* x, y compatible PADIC, op = add or sub */
616 static GEN
addsub_pp(GEN x,GEN y,GEN (* op)(GEN,GEN))617 addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
618 {
619 pari_sp av = avma;
620 long d,e,r,rx,ry;
621 GEN u, z, mod, p = gel(x,2);
622 int swap;
623
624 (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
625 e = valp(x);
626 r = valp(y); d = r-e;
627 if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
628 rx = precp(x);
629 ry = precp(y);
630 if (d) /* v(x) < v(y) */
631 {
632 r = d+ry; z = powiu(p,d);
633 if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
634 z = mulii(z,gel(y,4));
635 u = swap? op(z, gel(x,4)): op(gel(x,4), z);
636 }
637 else
638 {
639 long c;
640 if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
641 u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
642 if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
643 {
644 set_avma(av); return zeropadic(p, e+r);
645 }
646 if (c)
647 {
648 mod = diviiexact(mod, powiu(p,c));
649 r -= c;
650 e += c;
651 }
652 }
653 u = modii(u, mod);
654 set_avma(av); z = cgetg(5,t_PADIC);
655 z[1] = evalprecp(r) | evalvalp(e);
656 gel(z,2) = icopy(p);
657 gel(z,3) = icopy(mod);
658 gel(z,4) = icopy(u); return z;
659 }
660 /* Rg_to_Fp(t_FRAC) without GC */
661 static GEN
Q_to_Fp(GEN x,GEN p)662 Q_to_Fp(GEN x, GEN p)
663 { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
664 /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
665 static GEN
addQp(GEN x,GEN y)666 addQp(GEN x, GEN y)
667 {
668 pari_sp av = avma;
669 long d, r, e, vy = valp(y), py = precp(y);
670 GEN z, q, mod, u, p = gel(y,2);
671
672 e = Q_pvalrem(x, p, &x);
673 d = vy - e; r = d + py;
674 if (r <= 0) { set_avma(av); return gcopy(y); }
675 mod = gel(y,3);
676 u = gel(y,4);
677 (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
678
679 if (d > 0)
680 {
681 q = powiu(p,d);
682 mod = mulii(mod, q);
683 if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
684 u = addii(x, mulii(u, q));
685 }
686 else if (d < 0)
687 {
688 q = powiu(p,-d);
689 if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
690 u = addii(u, mulii(x, q));
691 r = py; e = vy;
692 }
693 else
694 {
695 long c;
696 if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
697 u = addii(u, x);
698 if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
699 {
700 set_avma(av); return zeropadic(p,e+r);
701 }
702 if (c)
703 {
704 mod = diviiexact(mod, powiu(p,c));
705 r -= c;
706 e += c;
707 }
708 }
709 u = modii(u, mod); set_avma(av);
710 z = cgetg(5,t_PADIC);
711 z[1] = evalprecp(r) | evalvalp(e);
712 gel(z,2) = icopy(p);
713 gel(z,3) = icopy(mod);
714 gel(z,4) = icopy(u); return z;
715 }
716
717 /* Mod(x,X) + Mod(y,X) */
718 #define addsub_polmod_same addsub_polmod_scal
719 /* Mod(x,X) +/- Mod(y,Y) */
720 static GEN
addsub_polmod(GEN X,GEN Y,GEN x,GEN y,GEN (* op)(GEN,GEN))721 addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
722 {
723 long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
724 GEN z = cgetg(3,t_POLMOD);
725 long vx = varn(X), vy = varn(Y);
726 if (vx==vy) {
727 pari_sp av;
728 gel(z,1) = RgX_gcd(X,Y); av = avma;
729 warn_coercion(X,Y,gel(z,1));
730 gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
731 }
732 if (varncmp(vx, vy) < 0)
733 { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
734 else
735 { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
736 gel(z,2) = op(x, y); return z;
737 }
738 /* Mod(y, Y) +/- x, x scalar or polynomial in same var and reduced degree */
739 static GEN
addsub_polmod_scal(GEN Y,GEN y,GEN x,GEN (* op)(GEN,GEN))740 addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
741 { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
742
743 /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
744 static GEN
add_ser_scal(GEN y,GEN x)745 add_ser_scal(GEN y, GEN x)
746 {
747 long i, l, ly, vy;
748 GEN z;
749
750 if (isrationalzero(x)) return gcopy(y);
751 ly = lg(y);
752 l = valp(y);
753 if (l < 3-ly) return gcopy(y);
754 /* l + ly >= 3 */
755 if (l < 0)
756 {
757 z = cgetg(ly,t_SER); z[1] = y[1];
758 for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
759 gel(z,i) = gadd(x,gel(y,i)); i++;
760 for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
761 return z;
762 }
763 vy = varn(y);
764 if (l > 0)
765 {
766 if (ser_isexactzero(y))
767 return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, l);
768 y -= l; ly += l;
769 z = cgetg(ly,t_SER);
770 x = gcopy(x);
771 for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
772 }
773 else
774 { /* l = 0, ly >= 3. Also OK if ser_isexactzero(y) */
775 z = cgetg(ly,t_SER);
776 x = gadd(x, gel(y,2));
777 i = 3;
778 }
779 for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
780 gel(z,2) = x;
781 z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
782 return gequal0(x)? normalize(z): z;
783 }
784 static long
_serprec(GEN x)785 _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
786 /* x,y t_SER in the same variable: x+y */
787 static GEN
ser_add(GEN x,GEN y)788 ser_add(GEN x, GEN y)
789 {
790 long i, lx,ly, n = valp(y) - valp(x);
791 GEN z;
792 if (n < 0) { n = -n; swap(x,y); }
793 /* valp(x) <= valp(y) */
794 lx = _serprec(x);
795 if (lx == 2) /* don't lose type information */
796 return scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), valp(x));
797 ly = _serprec(y) + n; if (lx < ly) ly = lx;
798 if (n)
799 {
800 if (n+2 > lx) return gcopy(x);
801 z = cgetg(ly,t_SER);
802 for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
803 for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
804 } else {
805 z = cgetg(ly,t_SER);
806 for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
807 }
808 z[1] = x[1]; return normalize(z);
809 }
810 /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
811 static GEN
add_rfrac_scal(GEN y,GEN x)812 add_rfrac_scal(GEN y, GEN x)
813 {
814 pari_sp av;
815 GEN n;
816
817 if (isintzero(x)) return gcopy(y); /* frequent special case */
818 av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
819 return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
820 }
821
822 /* x "scalar", ty != t_MAT and nonscalar */
823 static GEN
add_scal(GEN y,GEN x,long ty)824 add_scal(GEN y, GEN x, long ty)
825 {
826 switch(ty)
827 {
828 case t_POL: return RgX_Rg_add(y, x);
829 case t_SER: return add_ser_scal(y, x);
830 case t_RFRAC: return add_rfrac_scal(y, x);
831 case t_COL: return RgC_Rg_add(y, x);
832 case t_VEC:
833 if (isintzero(x)) return gcopy(y);
834 break;
835 }
836 pari_err_TYPE2("+",x,y);
837 return NULL; /* LCOV_EXCL_LINE */
838 }
839
840 /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
841 static GEN
setfrac(GEN z,GEN a,GEN b)842 setfrac(GEN z, GEN a, GEN b)
843 {
844 gel(z,1) = icopy_avma(a, (pari_sp)z);
845 gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
846 set_avma((pari_sp)gel(z,2)); return z;
847 }
848 /* z <- a / (b*Q), (Q,a) = 1 */
849 static GEN
addsub_frac_i(GEN z,GEN Q,GEN a,GEN b)850 addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
851 {
852 GEN q = Qdivii(a, b); /* != 0 */
853 if (typ(q) == t_INT)
854 {
855 gel(z,1) = gerepileuptoint((pari_sp)Q, q);
856 gel(z,2) = Q; return z;
857 }
858 return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
859 }
860 static GEN
addsub_frac(GEN x,GEN y,GEN (* op)(GEN,GEN))861 addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
862 {
863 GEN x1 = gel(x,1), x2 = gel(x,2);
864 GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
865 int s = cmpii(x2, y2);
866
867 /* common denominator: (x1 op y1) / x2 */
868 if (!s)
869 {
870 pari_sp av = avma;
871 return gerepileupto(av, Qdivii(op(x1, y1), x2));
872 }
873 z = cgetg(3, t_FRAC);
874 if (s < 0)
875 {
876 Q = dvmdii(y2, x2, &r);
877 /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
878 if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
879 delta = gcdii(x2,r);
880 }
881 else
882 {
883 Q = dvmdii(x2, y2, &r);
884 /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
885 if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
886 delta = gcdii(y2,r);
887 }
888 /* delta = gcd(x2,y2) */
889 if (equali1(delta))
890 { /* numerator is nonzero */
891 gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
892 gel(z,2) = mulii(x2,y2); return z;
893 }
894 x2 = diviiexact(x2,delta);
895 y2 = diviiexact(y2,delta);
896 n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
897 q = dvmdii(n, delta, &r);
898 if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
899 r = gcdii(delta, r);
900 if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
901 return setfrac(z, n, mulii(mulii(x2, y2), delta));
902 }
903
904 /* assume x2, y2 are t_POLs in the same variable */
905 static GEN
add_rfrac(GEN x,GEN y)906 add_rfrac(GEN x, GEN y)
907 {
908 pari_sp av = avma;
909 GEN x1 = gel(x,1), x2 = gel(x,2);
910 GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
911
912 delta = RgX_gcd(x2,y2);
913 if (!degpol(delta))
914 {
915 n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
916 d = RgX_mul(x2, y2);
917 return gerepileupto(av, gred_rfrac_simple(n, d));
918 }
919 x2 = RgX_div(x2,delta);
920 y2 = RgX_div(y2,delta);
921 n = gadd(gmul(x1,y2), gmul(y1,x2));
922 if (!signe(n))
923 {
924 n = simplify_shallow(n);
925 if (isrationalzero(n)) return gerepileupto(av, zeropol(varn(x2)));
926 return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
927 }
928 if (degpol(n) == 0)
929 return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
930 q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
931 if (isexactzero(r))
932 {
933 GEN z;
934 d = RgX_mul(x2, y2);
935 /* "constant" denominator ? */
936 z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
937 return gerepileupto(av, z);
938 }
939 r = RgX_gcd(delta, r);
940 if (degpol(r))
941 {
942 n = RgX_div(n, r);
943 d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
944 }
945 else
946 d = RgX_mul(gel(x,2), y2);
947 return gerepileupto(av, gred_rfrac_simple(n, d));
948 }
949
950 GEN
gadd(GEN x,GEN y)951 gadd(GEN x, GEN y)
952 {
953 long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
954 pari_sp av, tetpil;
955 GEN z, p1;
956
957 if (tx == ty) switch(tx) /* shortcut to generic case */
958 {
959 case t_INT: return addii(x,y);
960 case t_REAL: return addrr(x,y);
961 case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
962 z = cgetg(3,t_INTMOD);
963 if (X==Y || equalii(X,Y))
964 return add_intmod_same(z, X, gel(x,2), gel(y,2));
965 gel(z,1) = gcdii(X,Y);
966 warn_coercion(X,Y,gel(z,1));
967 av = avma; p1 = addii(gel(x,2),gel(y,2));
968 gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
969 }
970 case t_FRAC: return addsub_frac(x,y,addii);
971 case t_COMPLEX: z = cgetg(3,t_COMPLEX);
972 gel(z,2) = gadd(gel(x,2),gel(y,2));
973 if (isintzero(gel(z,2)))
974 {
975 set_avma((pari_sp)(z+3));
976 return gadd(gel(x,1),gel(y,1));
977 }
978 gel(z,1) = gadd(gel(x,1),gel(y,1));
979 return z;
980 case t_PADIC:
981 if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
982 return addsub_pp(x,y, addii);
983 case t_QUAD: z = cgetg(4,t_QUAD);
984 if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
985 gel(z,1) = ZX_copy(gel(x,1));
986 gel(z,2) = gadd(gel(x,2),gel(y,2));
987 gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
988 case t_POLMOD:
989 if (RgX_equal_var(gel(x,1), gel(y,1)))
990 return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
991 return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
992 case t_FFELT: return FF_add(x,y);
993 case t_POL:
994 vx = varn(x);
995 vy = varn(y);
996 if (vx != vy) {
997 if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
998 else return RgX_Rg_add(y, x);
999 }
1000 return RgX_add(x, y);
1001 case t_SER:
1002 vx = varn(x);
1003 vy = varn(y);
1004 if (vx != vy) {
1005 if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
1006 else return add_ser_scal(y, x);
1007 }
1008 return ser_add(x, y);
1009 case t_RFRAC:
1010 vx = varn(gel(x,2));
1011 vy = varn(gel(y,2));
1012 if (vx != vy) {
1013 if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
1014 else return add_rfrac_scal(y, x);
1015 }
1016 return add_rfrac(x,y);
1017 case t_VEC:
1018 if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1019 return RgV_add(x,y);
1020 case t_COL:
1021 if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1022 return RgC_add(x,y);
1023 case t_MAT:
1024 lx = lg(x);
1025 if (lg(y) != lx) pari_err_OP("+",x,y);
1026 if (lx == 1) return cgetg(1, t_MAT);
1027 if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1028 return RgM_add(x,y);
1029
1030 default: pari_err_TYPE2("+",x,y);
1031 }
1032 /* tx != ty */
1033 if (tx > ty) { swap(x,y); lswap(tx,ty); }
1034
1035 if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
1036 {
1037 case t_INT:
1038 switch(ty)
1039 {
1040 case t_REAL: return addir(x,y);
1041 case t_INTMOD:
1042 z = cgetg(3, t_INTMOD);
1043 return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1044 case t_FRAC: z = cgetg(3,t_FRAC);
1045 gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
1046 gel(z,2) = icopy(gel(y,2)); return z;
1047 case t_COMPLEX: return addRc(x, y);
1048 case t_PADIC:
1049 if (!signe(x)) return gcopy(y);
1050 return addQp(x,y);
1051 case t_QUAD: return addRq(x, y);
1052 case t_FFELT: return FF_Z_add(y,x);
1053 }
1054
1055 case t_REAL:
1056 switch(ty)
1057 {
1058 case t_FRAC:
1059 if (!signe(gel(y,1))) return rcopy(x);
1060 if (!signe(x))
1061 {
1062 lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
1063 return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
1064 }
1065 av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
1066 return gerepile(av,tetpil,divri(z,gel(y,2)));
1067 case t_COMPLEX: return addRc(x, y);
1068 case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
1069
1070 default: pari_err_TYPE2("+",x,y);
1071 }
1072
1073 case t_INTMOD:
1074 switch(ty)
1075 {
1076 case t_FRAC: { GEN X = gel(x,1);
1077 z = cgetg(3, t_INTMOD);
1078 p1 = Fp_div(gel(y,1), gel(y,2), X);
1079 return add_intmod_same(z, X, p1, gel(x,2));
1080 }
1081 case t_FFELT:
1082 if (!equalii(gel(x,1),FF_p_i(y)))
1083 pari_err_OP("+",x,y);
1084 return FF_Z_add(y,gel(x,2));
1085 case t_COMPLEX: return addRc(x, y);
1086 case t_PADIC: { GEN X = gel(x,1);
1087 z = cgetg(3, t_INTMOD);
1088 return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1089 }
1090 case t_QUAD: return addRq(x, y);
1091 }
1092
1093 case t_FRAC:
1094 switch (ty)
1095 {
1096 case t_COMPLEX: return addRc(x, y);
1097 case t_PADIC:
1098 if (!signe(gel(x,1))) return gcopy(y);
1099 return addQp(x,y);
1100 case t_QUAD: return addRq(x, y);
1101 case t_FFELT: return FF_Q_add(y, x);
1102 }
1103
1104 case t_FFELT:
1105 pari_err_TYPE2("+",x,y);
1106
1107 case t_COMPLEX:
1108 switch(ty)
1109 {
1110 case t_PADIC:
1111 return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
1112 case t_QUAD:
1113 lx = precision(x); if (!lx) pari_err_OP("+",x,y);
1114 return gequal0(y)? gcopy(x): addqf(y, x, lx);
1115 }
1116
1117 case t_PADIC: /* ty == t_QUAD */
1118 return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
1119 }
1120 /* tx < ty, !is_const_t(y) */
1121 switch(ty)
1122 {
1123 case t_MAT:
1124 if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
1125 if (isrationalzero(x)) return gcopy(y);
1126 return RgM_Rg_add(y, x);
1127 case t_COL:
1128 if (tx == t_VEC) pari_err_TYPE2("+",x,y);
1129 return RgC_Rg_add(y, x);
1130 case t_POLMOD: /* is_const_t(tx) in this case */
1131 return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
1132 }
1133 if (is_scalar_t(tx)) {
1134 if (tx == t_POLMOD)
1135 {
1136 vx = varn(gel(x,1));
1137 vy = gvar(y);
1138 if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
1139 else
1140 if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
1141 return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
1142 }
1143 return add_scal(y, x, ty);
1144 }
1145 /* x and y are not scalars, ty != t_MAT */
1146 vx = gvar(x);
1147 vy = gvar(y);
1148 if (vx != vy) { /* x or y is treated as a scalar */
1149 if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
1150 return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
1151 : add_scal(y, x, ty);
1152 }
1153 /* vx = vy */
1154 switch(tx)
1155 {
1156 case t_POL:
1157 switch (ty)
1158 {
1159 case t_SER:
1160 if (lg(x) == 2) return gcopy(y);
1161 i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
1162 i = lg(y) + valp(y) - i;
1163 if (i < 3) return gcopy(y);
1164 p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
1165 settyp(p1, t_VECSMALL); /* p1 left on stack */
1166 return y;
1167
1168 case t_RFRAC: return add_rfrac_scal(y, x);
1169 }
1170 break;
1171
1172 case t_SER:
1173 if (ty == t_RFRAC)
1174 {
1175 GEN n, d;
1176 long vn, vd;
1177 av = avma;
1178 n = gel(y,1); vn = gval(n, vy);
1179 d = gel(y,2); vd = RgX_valrem(d, &d);
1180
1181 l = lg(x) + valp(x) - (vn - vd);
1182 if (l < 3) { set_avma(av); return gcopy(x); }
1183
1184 /* take advantage of y = t^n ! */
1185 if (degpol(d))
1186 y = gdiv(n, RgX_to_ser_inexact(d,l));
1187 else {
1188 y = gdiv(n, gel(d,2));
1189 if (gvar(y) == vy) y = RgX_to_ser(y,l); else y = scalarser(y, vy, l);
1190 }
1191 setvalp(y, valp(y) - vd);
1192 return gerepileupto(av, gadd(y, x));
1193 }
1194 break;
1195 }
1196 pari_err_TYPE2("+",x,y);
1197 return NULL; /* LCOV_EXCL_LINE */
1198 }
1199
1200 GEN
gaddsg(long x,GEN y)1201 gaddsg(long x, GEN y)
1202 {
1203 long ty = typ(y);
1204 GEN z;
1205
1206 switch(ty)
1207 {
1208 case t_INT: return addsi(x,y);
1209 case t_REAL: return addsr(x,y);
1210 case t_INTMOD:
1211 z = cgetg(3, t_INTMOD);
1212 return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
1213 case t_FRAC: z = cgetg(3,t_FRAC);
1214 gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
1215 gel(z,2) = icopy(gel(y,2)); return z;
1216 case t_COMPLEX:
1217 z = cgetg(3, t_COMPLEX);
1218 gel(z,1) = gaddsg(x, gel(y,1));
1219 gel(z,2) = gcopy(gel(y,2)); return z;
1220
1221 default: return gadd(stoi(x), y);
1222 }
1223 }
1224
1225 GEN
gsubsg(long x,GEN y)1226 gsubsg(long x, GEN y)
1227 {
1228 GEN z, a, b;
1229 pari_sp av;
1230
1231 switch(typ(y))
1232 {
1233 case t_INT: return subsi(x,y);
1234 case t_REAL: return subsr(x,y);
1235 case t_INTMOD:
1236 z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
1237 return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
1238 case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
1239 gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
1240 gel(z,2) = icopy(gel(y,2)); return z;
1241 case t_COMPLEX:
1242 z = cgetg(3, t_COMPLEX);
1243 gel(z,1) = gsubsg(x, gel(y,1));
1244 gel(z,2) = gneg(gel(y,2)); return z;
1245 }
1246 av = avma;
1247 return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
1248 }
1249
1250 /********************************************************************/
1251 /** **/
1252 /** SUBTRACTION **/
1253 /** **/
1254 /********************************************************************/
1255
1256 GEN
gsub(GEN x,GEN y)1257 gsub(GEN x, GEN y)
1258 {
1259 long tx = typ(x), ty = typ(y);
1260 pari_sp av;
1261 GEN z;
1262 if (tx == ty) switch(tx) /* shortcut to generic case */
1263 {
1264 case t_INT: return subii(x,y);
1265 case t_REAL: return subrr(x,y);
1266 case t_INTMOD: { GEN p1, X = gel(x,1), Y = gel(y,1);
1267 z = cgetg(3,t_INTMOD);
1268 if (X==Y || equalii(X,Y))
1269 return sub_intmod_same(z, X, gel(x,2), gel(y,2));
1270 gel(z,1) = gcdii(X,Y);
1271 warn_coercion(X,Y,gel(z,1));
1272 av = avma; p1 = subii(gel(x,2),gel(y,2));
1273 gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
1274 }
1275 case t_FRAC: return addsub_frac(x,y, subii);
1276 case t_COMPLEX: z = cgetg(3,t_COMPLEX);
1277 gel(z,2) = gsub(gel(x,2),gel(y,2));
1278 if (isintzero(gel(z,2)))
1279 {
1280 set_avma((pari_sp)(z+3));
1281 return gsub(gel(x,1),gel(y,1));
1282 }
1283 gel(z,1) = gsub(gel(x,1),gel(y,1));
1284 return z;
1285 case t_PADIC:
1286 if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
1287 return addsub_pp(x,y, subii);
1288 case t_QUAD: z = cgetg(4,t_QUAD);
1289 if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
1290 gel(z,1) = ZX_copy(gel(x,1));
1291 gel(z,2) = gsub(gel(x,2),gel(y,2));
1292 gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
1293 case t_POLMOD:
1294 if (RgX_equal_var(gel(x,1), gel(y,1)))
1295 return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
1296 return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
1297 case t_FFELT: return FF_sub(x,y);
1298 case t_POL: {
1299 long vx = varn(x);
1300 long vy = varn(y);
1301 if (vx != vy) {
1302 if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
1303 else return Rg_RgX_sub(x, y);
1304 }
1305 return RgX_sub(x, y);
1306 }
1307 case t_VEC:
1308 if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1309 return RgV_sub(x,y);
1310 case t_COL:
1311 if (lg(y) != lg(x)) pari_err_OP("+",x,y);
1312 return RgC_sub(x,y);
1313 case t_MAT: {
1314 long lx = lg(x);
1315 if (lg(y) != lx) pari_err_OP("+",x,y);
1316 if (lx == 1) return cgetg(1, t_MAT);
1317 if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
1318 return RgM_sub(x,y);
1319 }
1320 case t_RFRAC: case t_SER: break;
1321
1322 default: pari_err_TYPE2("+",x,y);
1323 }
1324 av = avma;
1325 return gerepileupto(av, gadd(x,gneg_i(y)));
1326 }
1327
1328 /********************************************************************/
1329 /** **/
1330 /** MULTIPLICATION **/
1331 /** **/
1332 /********************************************************************/
1333 static GEN
mul_ser_scal(GEN y,GEN x)1334 mul_ser_scal(GEN y, GEN x) {
1335 long l, i;
1336 GEN z;
1337 if (isrationalzero(x)) return gmul(Rg_get_0(y), x);
1338 if (ser_isexactzero(y))
1339 return scalarser(lg(y) == 2? Rg_get_0(x)
1340 : gmul(gel(y,2), x), varn(y), valp(y));
1341 z = cgetg_copy(y, &l); z[1] = y[1];
1342 for (i = 2; i < l; i++) gel(z,i) = gmul(gel(y,i), x);
1343 return normalize(z);
1344 }
1345 /* (n/d) * x, x "scalar" or polynomial in the same variable as d
1346 * [n/d a valid RFRAC] */
1347 static GEN
mul_rfrac_scal(GEN n,GEN d,GEN x)1348 mul_rfrac_scal(GEN n, GEN d, GEN x)
1349 {
1350 pari_sp av = avma;
1351 GEN z;
1352
1353 switch(typ(x))
1354 {
1355 case t_PADIC:
1356 n = gmul(n, x);
1357 d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
1358 return gerepileupto(av, gdiv(n,d));
1359
1360 case t_INTMOD: case t_POLMOD:
1361 n = gmul(n, x);
1362 d = gmul(d, gmodulo(gen_1, gel(x,1)));
1363 return gerepileupto(av, gdiv(n,d));
1364 }
1365 z = gred_rfrac2(x, d);
1366 n = simplify_shallow(n);
1367 if (typ(z) == t_RFRAC)
1368 {
1369 n = gmul(gel(z,1), n);
1370 d = gel(z,2);
1371 if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
1372 z = RgX_Rg_div(n, d);
1373 else
1374 z = gred_rfrac_simple(n, d);
1375 }
1376 else
1377 z = gmul(z, n);
1378 return gerepileupto(av, z);
1379 }
1380 static GEN
mul_scal(GEN y,GEN x,long ty)1381 mul_scal(GEN y, GEN x, long ty)
1382 {
1383 switch(ty)
1384 {
1385 case t_POL:
1386 if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
1387 return RgX_Rg_mul(y, x);
1388 case t_SER: return mul_ser_scal(y, x);
1389 case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1390 case t_QFI: case t_QFR:
1391 if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
1392 }
1393 pari_err_TYPE2("*",x,y);
1394 return NULL; /* LCOV_EXCL_LINE */
1395 }
1396
1397 static GEN
mul_gen_rfrac(GEN X,GEN Y)1398 mul_gen_rfrac(GEN X, GEN Y)
1399 {
1400 GEN y1 = gel(Y,1), y2 = gel(Y,2);
1401 long vx = gvar(X), vy = varn(y2);
1402 return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1403 gred_rfrac_simple(gmul(y1,X), y2);
1404 }
1405 /* (x1/x2) * (y1/y2) */
1406 static GEN
mul_rfrac(GEN x1,GEN x2,GEN y1,GEN y2)1407 mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1408 {
1409 GEN z, X, Y;
1410 pari_sp av = avma;
1411
1412 X = gred_rfrac2(x1, y2);
1413 Y = gred_rfrac2(y1, x2);
1414 if (typ(X) == t_RFRAC)
1415 {
1416 if (typ(Y) == t_RFRAC) {
1417 x1 = gel(X,1);
1418 x2 = gel(X,2);
1419 y1 = gel(Y,1);
1420 y2 = gel(Y,2);
1421 z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1422 } else
1423 z = mul_gen_rfrac(Y, X);
1424 }
1425 else if (typ(Y) == t_RFRAC)
1426 z = mul_gen_rfrac(X, Y);
1427 else
1428 z = gmul(X, Y);
1429 return gerepileupto(av, z);
1430 }
1431 /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
1432 static GEN
div_rfrac_pol(GEN x1,GEN x2,GEN y2)1433 div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1434 {
1435 pari_sp av = avma;
1436 GEN X = gred_rfrac2(x1, y2);
1437 if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
1438 {
1439 x2 = RgX_mul(gel(X,2), x2);
1440 x1 = gel(X,1);
1441 }
1442 else
1443 x1 = X;
1444 return gerepileupto(av, gred_rfrac_simple(x1, x2));
1445 }
1446
1447 /* Mod(y, Y) * x, assuming x scalar */
1448 static GEN
mul_polmod_scal(GEN Y,GEN y,GEN x)1449 mul_polmod_scal(GEN Y, GEN y, GEN x)
1450 { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
1451
1452 /* cf mulqq */
1453 static GEN
quad_polmod_mul(GEN P,GEN x,GEN y)1454 quad_polmod_mul(GEN P, GEN x, GEN y)
1455 {
1456 GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
1457 pari_sp tetpil, av = avma;
1458 T[1] = x[1];
1459 p2 = gmul(gel(x,2), gel(y,2));
1460 p3 = gmul(gel(x,3), gel(y,3));
1461 p1 = gmul(gneg_i(c),p3);
1462 /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
1463 if (typ(b) == t_INT)
1464 {
1465 if (signe(b))
1466 {
1467 p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1468 if (is_pm1(b))
1469 {
1470 if (signe(b) > 0) p3 = gneg(p3);
1471 }
1472 else
1473 p3 = gmul(negi(b), p3);
1474 }
1475 else
1476 {
1477 p3 = gmul(gel(x,2),gel(y,3));
1478 p4 = gmul(gel(x,3),gel(y,2));
1479 }
1480 }
1481 else
1482 {
1483 p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
1484 p3 = gmul(gneg_i(b), p3);
1485 }
1486 tetpil = avma;
1487 gel(T,2) = gadd(p2, p1);
1488 gel(T,3) = gadd(p4, p3);
1489 gerepilecoeffssp(av,tetpil,T+2,2);
1490 return normalizepol_lg(T,4);
1491 }
1492 /* Mod(x,T) * Mod(y,T) */
1493 static GEN
mul_polmod_same(GEN T,GEN x,GEN y)1494 mul_polmod_same(GEN T, GEN x, GEN y)
1495 {
1496 GEN z = cgetg(3,t_POLMOD), a;
1497 long v = varn(T), lx = lg(x), ly = lg(y);
1498 gel(z,1) = RgX_copy(T);
1499 /* x * y mod T optimised */
1500 if (typ(x) != t_POL || varn(x) != v || lx <= 3
1501 || typ(y) != t_POL || varn(y) != v || ly <= 3)
1502 a = gmul(x, y);
1503 else
1504 {
1505 if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
1506 a = quad_polmod_mul(T, x, y);
1507 else
1508 a = RgXQ_mul(x, y, gel(z,1));
1509 }
1510 gel(z,2) = a; return z;
1511 }
1512 static GEN
sqr_polmod(GEN T,GEN x)1513 sqr_polmod(GEN T, GEN x)
1514 {
1515 GEN a, z = cgetg(3,t_POLMOD);
1516 gel(z,1) = RgX_copy(T);
1517 if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
1518 a = gsqr(x);
1519 else
1520 {
1521 pari_sp av = avma;
1522 a = RgXQ_sqr(x, gel(z,1));
1523 a = gerepileupto(av, a);
1524 }
1525 gel(z,2) = a; return z;
1526 }
1527 /* Mod(x,X) * Mod(y,Y) */
1528 static GEN
mul_polmod(GEN X,GEN Y,GEN x,GEN y)1529 mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1530 {
1531 long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1532 long vx = varn(X), vy = varn(Y);
1533 GEN z = cgetg(3,t_POLMOD);
1534
1535 if (vx==vy) {
1536 pari_sp av;
1537 gel(z,1) = RgX_gcd(X,Y); av = avma;
1538 warn_coercion(X,Y,gel(z,1));
1539 gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
1540 return z;
1541 }
1542 if (varncmp(vx, vy) < 0)
1543 { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1544 else
1545 { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1546 gel(z,2) = gmul(x, y); return z;
1547 }
1548
1549 #if 0 /* used by 3M only */
1550 /* set z = x+y and return 1 if x,y have the same sign
1551 * set z = x-y and return 0 otherwise */
1552 static int
1553 did_add(GEN x, GEN y, GEN *z)
1554 {
1555 long tx = typ(x), ty = typ(y);
1556 if (tx == ty) switch(tx)
1557 {
1558 case t_INT: *z = addii(x,y); return 1;
1559 case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
1560 case t_REAL:
1561 if (signe(x) == -signe(y))
1562 { *z = subrr(x,y); return 0; }
1563 else
1564 { *z = addrr(x,y); return 1; }
1565 }
1566 if (tx == t_REAL) switch(ty)
1567 {
1568 case t_INT:
1569 if (signe(x) == -signe(y))
1570 { *z = subri(x,y); return 0; }
1571 else
1572 { *z = addri(x,y); return 1; }
1573 case t_FRAC:
1574 if (signe(x) == -signe(gel(y,1)))
1575 { *z = gsub(x,y); return 0; }
1576 else
1577 { *z = gadd(x,y); return 1; }
1578 }
1579 else if (ty == t_REAL) switch(tx)
1580 {
1581 case t_INT:
1582 if (signe(x) == -signe(y))
1583 { *z = subir(x,y); return 0; }
1584 else
1585 { *z = addir(x,y); return 1; }
1586 case t_FRAC:
1587 if (signe(gel(x,1)) == -signe(y))
1588 { *z = gsub(x,y); return 0; }
1589 else
1590 { *z = gadd(x,y); return 1; }
1591 }
1592 *z = gadd(x,y); return 1;
1593 }
1594 #endif
1595 /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
1596 static GEN
mulcIR(GEN x,GEN y)1597 mulcIR(GEN x, GEN y)
1598 {
1599 GEN z = cgetg(3,t_COMPLEX);
1600 pari_sp av = avma;
1601 gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
1602 gel(z,2) = gmul(y, gel(x,1));
1603 return z;
1604
1605 }
1606 /* x,y COMPLEX */
1607 static GEN
mulcc(GEN x,GEN y)1608 mulcc(GEN x, GEN y)
1609 {
1610 GEN xr = gel(x,1), xi = gel(x,2);
1611 GEN yr = gel(y,1), yi = gel(y,2);
1612 GEN p1, p2, p3, p4, z;
1613 pari_sp tetpil, av;
1614
1615 if (isintzero(xr))
1616 {
1617 if (isintzero(yr)) {
1618 av = avma;
1619 return gerepileupto(av, gneg(gmul(xi,yi)));
1620 }
1621 return mulcIR(y, xi);
1622 }
1623 if (isintzero(yr)) return mulcIR(x, yi);
1624
1625 z = cgetg(3,t_COMPLEX); av = avma;
1626 #if 0
1627 /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
1628 * e.g. xr + xi if exponents differ */
1629 if (did_add(xr, xi, &p3))
1630 {
1631 if (did_add(yr, yi, &p4)) {
1632 /* R = xr*yr - xi*yi
1633 * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
1634 p1 = gmul(xr,yr);
1635 p2 = gmul(xi,yi); p2 = gneg(p2);
1636 p3 = gmul(p3, p4);
1637 p4 = gsub(p2, p1);
1638 } else {
1639 /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
1640 * I = xr*yi + xi*yr */
1641 p1 = gmul(p3,p4);
1642 p3 = gmul(xr,yi);
1643 p4 = gmul(xi,yr);
1644 p2 = gsub(p3, p4);
1645 }
1646 } else {
1647 if (did_add(yr, yi, &p4)) {
1648 /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
1649 * I = xr*yi +xi*yr */
1650 p1 = gmul(p3,p4);
1651 p3 = gmul(xr,yi);
1652 p4 = gmul(xi,yr);
1653 p2 = gsub(p4, p3);
1654 } else {
1655 /* R = xr*yr - xi*yi
1656 * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
1657 p3 = gneg( gmul(p3, p4) );
1658 p1 = gmul(xr,yr);
1659 p2 = gmul(xi,yi);
1660 p4 = gadd(p1, p2);
1661
1662 p2 = gneg(p2);
1663 }
1664 }
1665 tetpil = avma;
1666 gel(z,1) = gadd(p1,p2);
1667 gel(z,2) = gadd(p3,p4);
1668 #else
1669 if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
1670 { /* 3M formula */
1671 p3 = addii(xr,xi);
1672 p4 = addii(yr,yi);
1673 p1 = mulii(xr,yr);
1674 p2 = mulii(xi,yi);
1675 p3 = mulii(p3,p4);
1676 p4 = addii(p2,p1);
1677 tetpil = avma;
1678 gel(z,1) = subii(p1,p2);
1679 gel(z,2) = subii(p3,p4);
1680 if (!signe(gel(z,2)))
1681 return gerepileuptoint((pari_sp)(z+3), gel(z,1));
1682 }
1683 else
1684 { /* naive 4M formula: avoid all loss of accuracy */
1685 p1 = gmul(xr,yr);
1686 p2 = gmul(xi,yi);
1687 p3 = gmul(xr,yi);
1688 p4 = gmul(xi,yr);
1689 tetpil = avma;
1690 gel(z,1) = gsub(p1,p2);
1691 gel(z,2) = gadd(p3,p4);
1692 if (isintzero(gel(z,2)))
1693 {
1694 cgiv(gel(z,2));
1695 return gerepileupto((pari_sp)(z+3), gel(z,1));
1696 }
1697 }
1698 #endif
1699
1700 gerepilecoeffssp(av,tetpil, z+1,2); return z;
1701 }
1702 /* x,y PADIC */
1703 static GEN
mulpp(GEN x,GEN y)1704 mulpp(GEN x, GEN y) {
1705 long l = valp(x) + valp(y);
1706 pari_sp av;
1707 GEN z, t;
1708 if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
1709 if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
1710 if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
1711
1712 t = (precp(x) > precp(y))? y: x;
1713 z = cgetp(t); setvalp(z,l); av = avma;
1714 affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
1715 set_avma(av); return z;
1716 }
1717 /* x,y QUAD */
1718 static GEN
mulqq(GEN x,GEN y)1719 mulqq(GEN x, GEN y) {
1720 GEN z = cgetg(4,t_QUAD);
1721 GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
1722 pari_sp av, tetpil;
1723 if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
1724
1725 gel(z,1) = ZX_copy(P); av = avma;
1726 p2 = gmul(gel(x,2),gel(y,2));
1727 p3 = gmul(gel(x,3),gel(y,3));
1728 p1 = gmul(gneg_i(c),p3);
1729
1730 if (signe(b))
1731 p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
1732 else
1733 {
1734 p3 = gmul(gel(x,2),gel(y,3));
1735 p4 = gmul(gel(x,3),gel(y,2));
1736 }
1737 tetpil = avma;
1738 gel(z,2) = gadd(p2,p1);
1739 gel(z,3) = gadd(p4,p3);
1740 gerepilecoeffssp(av,tetpil,z+2,2); return z;
1741 }
1742
1743 GEN
mulcxI(GEN x)1744 mulcxI(GEN x)
1745 {
1746 GEN z;
1747 switch(typ(x))
1748 {
1749 case t_INT: case t_REAL: case t_FRAC:
1750 return mkcomplex(gen_0, x);
1751 case t_COMPLEX:
1752 if (isintzero(gel(x,1))) return gneg(gel(x,2));
1753 z = cgetg(3,t_COMPLEX);
1754 gel(z,1) = gneg(gel(x,2));
1755 gel(z,2) = gel(x,1); return z;
1756 default:
1757 return gmul(gen_I(), x);
1758 }
1759 }
1760 GEN
mulcxmI(GEN x)1761 mulcxmI(GEN x)
1762 {
1763 GEN z;
1764 switch(typ(x))
1765 {
1766 case t_INT: case t_REAL: case t_FRAC:
1767 return mkcomplex(gen_0, gneg(x));
1768 case t_COMPLEX:
1769 if (isintzero(gel(x,1))) return gel(x,2);
1770 z = cgetg(3,t_COMPLEX);
1771 gel(z,1) = gel(x,2);
1772 gel(z,2) = gneg(gel(x,1)); return z;
1773 default:
1774 return gmul(mkcomplex(gen_0, gen_m1), x);
1775 }
1776 }
1777 /* x * I^k */
1778 GEN
mulcxpowIs(GEN x,long k)1779 mulcxpowIs(GEN x, long k)
1780 {
1781 switch (k & 3)
1782 {
1783 case 1: return mulcxI(x);
1784 case 2: return gneg(x);
1785 case 3: return mulcxmI(x);
1786 }
1787 return x;
1788 }
1789
1790 /* fill in coefficients of t_SER z from coeffs of t_POL y */
1791 static GEN
fill_ser(GEN z,GEN y)1792 fill_ser(GEN z, GEN y)
1793 {
1794 long i, lx = lg(z), ly = lg(y);
1795 if (ly >= lx) {
1796 for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
1797 } else {
1798 for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
1799 for ( ; i < lx; i++) gel(z,i) = gen_0;
1800 }
1801 return normalize(z);
1802 }
1803
1804 GEN
gmul(GEN x,GEN y)1805 gmul(GEN x, GEN y)
1806 {
1807 long tx, ty, lx, ly, vx, vy, i, l;
1808 pari_sp av, tetpil;
1809 GEN z, p1, p2;
1810
1811 if (x == y) return gsqr(x);
1812 tx = typ(x); ty = typ(y);
1813 if (tx == ty) switch(tx)
1814 {
1815 case t_INT: return mulii(x,y);
1816 case t_REAL: return mulrr(x,y);
1817 case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1818 z = cgetg(3,t_INTMOD);
1819 if (X==Y || equalii(X,Y))
1820 return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1821 gel(z,1) = gcdii(X,Y);
1822 warn_coercion(X,Y,gel(z,1));
1823 av = avma; p1 = mulii(gel(x,2),gel(y,2));
1824 gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
1825 }
1826 case t_FRAC:
1827 {
1828 GEN x1 = gel(x,1), x2 = gel(x,2);
1829 GEN y1 = gel(y,1), y2 = gel(y,2);
1830 z=cgetg(3,t_FRAC);
1831 p1 = gcdii(x1, y2);
1832 if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1833 p1 = gcdii(x2, y1);
1834 if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1835 tetpil = avma;
1836 gel(z,2) = mulii(x2,y2);
1837 gel(z,1) = mulii(x1,y1);
1838 fix_frac_if_int_GC(z,tetpil); return z;
1839 }
1840 case t_COMPLEX: return mulcc(x, y);
1841 case t_PADIC: return mulpp(x, y);
1842 case t_QUAD: return mulqq(x, y);
1843 case t_FFELT: return FF_mul(x, y);
1844 case t_POLMOD:
1845 if (RgX_equal_var(gel(x,1), gel(y,1)))
1846 return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1847 return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1848 case t_POL:
1849 vx = varn(x);
1850 vy = varn(y);
1851 if (vx != vy) {
1852 if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1853 else return RgX_Rg_mul(y, x);
1854 }
1855 return RgX_mul(x, y);
1856
1857 case t_SER: {
1858 vx = varn(x);
1859 vy = varn(y);
1860 if (vx != vy) {
1861 if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1862 return mul_ser_scal(y, x);
1863 }
1864 lx = minss(lg(x), lg(y));
1865 if (lx == 2) return zeroser(vx, valp(x)+valp(y));
1866 av = avma; z = cgetg(lx,t_SER);
1867 z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
1868 x = ser2pol_i(x, lx);
1869 y = ser2pol_i(y, lx);
1870 y = RgXn_mul(x, y, lx-2);
1871 return gerepilecopy(av, fill_ser(z,y));
1872 }
1873 case t_QFI: return qficomp(x,y);
1874 case t_QFR: return qfrcomp(x,y);
1875 case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1876 case t_MAT: return RgM_mul(x, y);
1877
1878 case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1879 z = cgetg_copy(x, &l);
1880 if (l != lg(y)) break;
1881 for (i=1; i<l; i++)
1882 {
1883 long yi = y[i];
1884 if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
1885 z[i] = x[yi];
1886 }
1887 return z;
1888
1889 default:
1890 pari_err_TYPE2("*",x,y);
1891 }
1892 /* tx != ty */
1893 if (is_const_t(ty) && is_const_t(tx)) {
1894 if (tx > ty) { swap(x,y); lswap(tx,ty); }
1895 switch(tx) {
1896 case t_INT:
1897 switch(ty)
1898 {
1899 case t_REAL: return signe(x)? mulir(x,y): gen_0;
1900 case t_INTMOD:
1901 z = cgetg(3, t_INTMOD);
1902 return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1903 case t_FRAC:
1904 if (!signe(x)) return gen_0;
1905 z=cgetg(3,t_FRAC);
1906 p1 = gcdii(x,gel(y,2));
1907 if (equali1(p1))
1908 {
1909 set_avma((pari_sp)z);
1910 gel(z,2) = icopy(gel(y,2));
1911 gel(z,1) = mulii(gel(y,1), x);
1912 }
1913 else
1914 {
1915 x = diviiexact(x,p1); tetpil = avma;
1916 gel(z,2) = diviiexact(gel(y,2), p1);
1917 gel(z,1) = mulii(gel(y,1), x);
1918 fix_frac_if_int_GC(z,tetpil);
1919 }
1920 return z;
1921 case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
1922 case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1923 case t_QUAD: return mulRq(x,y);
1924 case t_FFELT: return FF_Z_mul(y,x);
1925 }
1926
1927 case t_REAL:
1928 switch(ty)
1929 {
1930 case t_FRAC: return mulrfrac(x, y);
1931 case t_COMPLEX: return mulRc(x, y);
1932 case t_QUAD: return mulqf(y, x, realprec(x));
1933 default: pari_err_TYPE2("*",x,y);
1934 }
1935
1936 case t_INTMOD:
1937 switch(ty)
1938 {
1939 case t_FRAC: { GEN X = gel(x,1);
1940 z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
1941 return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1942 }
1943 case t_COMPLEX: return mulRc_direct(x,y);
1944 case t_PADIC: { GEN X = gel(x,1);
1945 z = cgetg(3, t_INTMOD);
1946 return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1947 }
1948 case t_QUAD: return mulRq(x, y);
1949 case t_FFELT:
1950 if (!equalii(gel(x,1),FF_p_i(y)))
1951 pari_err_OP("*",x,y);
1952 return FF_Z_mul(y,gel(x,2));
1953 }
1954
1955 case t_FRAC:
1956 switch(ty)
1957 {
1958 case t_COMPLEX: return mulRc(x, y);
1959 case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
1960 case t_QUAD: return mulRq(x, y);
1961 case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
1962 }
1963
1964 case t_FFELT:
1965 pari_err_TYPE2("*",x,y);
1966
1967 case t_COMPLEX:
1968 switch(ty)
1969 {
1970 case t_PADIC:
1971 return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
1972 case t_QUAD:
1973 lx = precision(x); if (!lx) pari_err_OP("*",x,y);
1974 return mulqf(y, x, lx);
1975 }
1976
1977 case t_PADIC: /* ty == t_QUAD */
1978 return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
1979 }
1980 }
1981
1982 if (is_matvec_t(ty))
1983 {
1984 if (!is_matvec_t(tx))
1985 {
1986 if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
1987 z = cgetg_copy(y, &ly);
1988 for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
1989 return z;
1990 }
1991 switch(tx)
1992 {
1993 case t_VEC:
1994 switch(ty) {
1995 case t_COL: return RgV_RgC_mul(x,y);
1996 case t_MAT: return RgV_RgM_mul(x,y);
1997 }
1998 break;
1999 case t_COL:
2000 switch(ty) {
2001 case t_VEC: return RgC_RgV_mul(x,y);
2002 case t_MAT: return RgC_RgM_mul(x,y);
2003 }
2004 break;
2005 case t_MAT:
2006 switch(ty) {
2007 case t_VEC: return RgM_RgV_mul(x,y);
2008 case t_COL: return RgM_RgC_mul(x,y);
2009 }
2010 }
2011 }
2012 if (is_matvec_t(tx))
2013 {
2014 if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
2015 z = cgetg_copy(x, &lx);
2016 for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
2017 return z;
2018 }
2019 if (tx > ty) { swap(x,y); lswap(tx,ty); }
2020 /* tx < ty, !ismatvec(x and y) */
2021
2022 if (ty == t_POLMOD) /* is_const_t(tx) in this case */
2023 return mul_polmod_scal(gel(y,1), gel(y,2), x);
2024 if (is_scalar_t(tx)) {
2025 if (tx == t_POLMOD) {
2026 vx = varn(gel(x,1));
2027 vy = gvar(y);
2028 if (vx != vy) {
2029 if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
2030 return mul_polmod_scal(gel(x,1), gel(x,2), y);
2031 }
2032 /* error if ty == t_SER */
2033 av = avma; y = gmod(y, gel(x,1));
2034 return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
2035 }
2036 return mul_scal(y, x, ty);
2037 }
2038
2039 /* x and y are not scalars, nor matvec */
2040 vx = gvar(x);
2041 vy = gvar(y);
2042 if (vx != vy) /* x or y is treated as a scalar */
2043 return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
2044 : mul_scal(y, x, ty);
2045 /* vx = vy */
2046 switch(tx)
2047 {
2048 case t_POL:
2049 switch (ty)
2050 {
2051 case t_SER:
2052 {
2053 long vn;
2054 if (lg(x) == 2) return zeropol(vx);
2055 if (lg(y) == 2) return zeroser(vx, valp(y)+RgX_val(x));
2056 av = avma;
2057 vn = RgX_valrem(x, &x);
2058 /* take advantage of x = t^n ! */
2059 if (degpol(x)) {
2060 p1 = RgX_to_ser(x,lg(y));
2061 if (vn) settyp(x, t_VECSMALL); /* *new* x left on stack */
2062 p2 = gmul(p1,y);
2063 settyp(p1, t_VECSMALL); /* p1 left on stack */
2064 } else {
2065 set_avma(av);
2066 p2 = mul_ser_scal(y, gel(x,2));
2067 }
2068 setvalp(p2, valp(p2) + vn);
2069 return p2;
2070 }
2071
2072 case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
2073 }
2074 break;
2075
2076 case t_SER:
2077 switch (ty)
2078 {
2079 case t_RFRAC:
2080 av = avma;
2081 return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
2082 }
2083 break;
2084 }
2085 pari_err_TYPE2("*",x,y);
2086 return NULL; /* LCOV_EXCL_LINE */
2087 }
2088
2089 /* return a nonnormalized result */
2090 GEN
sqr_ser_part(GEN x,long l1,long l2)2091 sqr_ser_part(GEN x, long l1, long l2)
2092 {
2093 long i, j, l;
2094 pari_sp av;
2095 GEN Z, z, p1, p2;
2096 long mi;
2097 if (l2 < l1) return zeroser(varn(x), 2*valp(x));
2098 p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
2099 Z = cgetg(l2-l1+3, t_SER);
2100 Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
2101 z = Z + 2-l1;
2102 x += 2; mi = 0;
2103 for (i=0; i<l1; i++)
2104 {
2105 p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2106 }
2107
2108 for (i=l1; i<=l2; i++)
2109 {
2110 p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
2111 p1=gen_0; av=avma; l=((i+1)>>1) - 1;
2112 for (j=i-mi; j<=minss(l,mi); j++)
2113 if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
2114 p1 = gshift(p1,1);
2115 if ((i&1) == 0 && p2[i>>1])
2116 p1 = gadd(p1, gsqr(gel(x,i>>1)));
2117 gel(z,i) = gerepileupto(av,p1);
2118 }
2119 return Z;
2120 }
2121
2122 GEN
gsqr(GEN x)2123 gsqr(GEN x)
2124 {
2125 long i, lx;
2126 pari_sp av, tetpil;
2127 GEN z, p1, p2, p3, p4;
2128
2129 switch(typ(x))
2130 {
2131 case t_INT: return sqri(x);
2132 case t_REAL: return sqrr(x);
2133 case t_INTMOD: { GEN X = gel(x,1);
2134 z = cgetg(3,t_INTMOD);
2135 gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
2136 gel(z,1) = icopy(X); return z;
2137 }
2138 case t_FRAC: return sqrfrac(x);
2139
2140 case t_COMPLEX:
2141 if (isintzero(gel(x,1))) {
2142 av = avma;
2143 return gerepileupto(av, gneg(gsqr(gel(x,2))));
2144 }
2145 z = cgetg(3,t_COMPLEX); av = avma;
2146 p1 = gadd(gel(x,1),gel(x,2));
2147 p2 = gsub(gel(x,1), gel(x,2));
2148 p3 = gmul(gel(x,1),gel(x,2));
2149 tetpil = avma;
2150 gel(z,1) = gmul(p1,p2);
2151 gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
2152
2153 case t_PADIC:
2154 z = cgetg(5,t_PADIC);
2155 i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
2156 if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
2157 z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
2158 gel(z,2) = icopy(gel(x,2));
2159 gel(z,3) = shifti(gel(x,3), i); av = avma;
2160 gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
2161 return z;
2162
2163 case t_QUAD: z = cgetg(4,t_QUAD);
2164 p1 = gel(x,1);
2165 gel(z,1) = ZX_copy(p1); av = avma;
2166 p2 = gsqr(gel(x,2));
2167 p3 = gsqr(gel(x,3));
2168 p4 = gmul(gneg_i(gel(p1,2)),p3);
2169
2170 if (gequal0(gel(p1,3)))
2171 {
2172 tetpil = avma;
2173 gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
2174 av = avma;
2175 p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
2176 gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
2177 }
2178
2179 p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
2180 tetpil = avma;
2181 gel(z,2) = gadd(p2,p4);
2182 gel(z,3) = gadd(p1,p3);
2183 gerepilecoeffssp(av,tetpil,z+2,2); return z;
2184
2185 case t_POLMOD:
2186 return sqr_polmod(gel(x,1), gel(x,2));
2187
2188 case t_FFELT: return FF_sqr(x);
2189
2190 case t_POL: return RgX_sqr(x);
2191
2192 case t_SER:
2193 lx = lg(x);
2194 if (ser_isexactzero(x)) {
2195 GEN z = gcopy(x);
2196 setvalp(z, 2*valp(x));
2197 return z;
2198 }
2199 if (lx < 40)
2200 return normalize( sqr_ser_part(x, 0, lx-3) );
2201 else
2202 {
2203 pari_sp av = avma;
2204 GEN z = cgetg(lx,t_SER);
2205 z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x)) | evalsigne(1);
2206 x = ser2pol_i(x,lx);
2207 x = RgXn_sqr(x, lx-2);
2208 return gerepilecopy(av, fill_ser(z,x));
2209 }
2210
2211 case t_RFRAC: z = cgetg(3,t_RFRAC);
2212 gel(z,1) = gsqr(gel(x,1));
2213 gel(z,2) = gsqr(gel(x,2)); return z;
2214
2215 case t_MAT: return RgM_sqr(x);
2216 case t_QFR: return qfrsqr(x);
2217 case t_QFI: return qfisqr(x);
2218 case t_VECSMALL:
2219 z = cgetg_copy(x, &lx);
2220 for (i=1; i<lx; i++)
2221 {
2222 long xi = x[i];
2223 if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
2224 z[i] = x[xi];
2225 }
2226 return z;
2227 }
2228 pari_err_TYPE2("*",x,x);
2229 return NULL; /* LCOV_EXCL_LINE */
2230 }
2231
2232 /********************************************************************/
2233 /** **/
2234 /** DIVISION **/
2235 /** **/
2236 /********************************************************************/
2237 static GEN
div_rfrac_scal(GEN x,GEN y)2238 div_rfrac_scal(GEN x, GEN y)
2239 {
2240 pari_sp av = avma;
2241 GEN d = rfrac_denom_mul_scal(gel(x,2), y);
2242 return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
2243 }
2244 static GEN
div_scal_rfrac(GEN x,GEN y)2245 div_scal_rfrac(GEN x, GEN y)
2246 {
2247 GEN y1 = gel(y,1), y2 = gel(y,2);
2248 pari_sp av = avma;
2249 if (typ(y1) == t_POL && varn(y2) == varn(y1))
2250 {
2251 if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
2252 y1 = gel(y1,2);
2253 }
2254 return RgX_Rg_mul(y2, gdiv(x,y1));
2255 }
2256 static GEN
div_rfrac(GEN x,GEN y)2257 div_rfrac(GEN x, GEN y)
2258 { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
2259
2260 /* x != 0 */
2261 static GEN
div_ser_scal(GEN y,GEN x)2262 div_ser_scal(GEN y, GEN x) {
2263 long i, l;
2264 GEN z;
2265 if (ser_isexactzero(y))
2266 return scalarser(lg(y) == 2? Rg_get_0(x)
2267 : gdiv(gel(y,2), x), varn(y), valp(y));
2268 z = cgetg_copy(y, &l); z[1] = y[1];
2269 for (i = 2; i < l; i++) gel(z,i) = gdiv(gel(y,i), x);
2270 return normalize(z);
2271 }
2272 GEN
ser_normalize(GEN x)2273 ser_normalize(GEN x)
2274 {
2275 long i, lx = lg(x);
2276 GEN c, z;
2277 if (lx == 2) return x;
2278 c = gel(x,2); if (gequal1(c)) return x;
2279 z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
2280 for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
2281 return z;
2282 }
2283
2284 /* y != 0 */
2285 static GEN
div_T_scal(GEN x,GEN y,long tx)2286 div_T_scal(GEN x, GEN y, long tx) {
2287 switch(tx)
2288 {
2289 case t_POL: return RgX_Rg_div(x, y);
2290 case t_SER: return div_ser_scal(x, y);
2291 case t_RFRAC: return div_rfrac_scal(x,y);
2292 }
2293 pari_err_TYPE2("/",x,y);
2294 return NULL; /* LCOV_EXCL_LINE */
2295 }
2296
2297 static GEN
div_scal_pol(GEN x,GEN y)2298 div_scal_pol(GEN x, GEN y) {
2299 long ly = lg(y);
2300 pari_sp av;
2301 if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
2302 if (isrationalzero(x)) return zeropol(varn(y));
2303 av = avma;
2304 return gerepileupto(av, gred_rfrac_simple(x,y));
2305 }
2306 static GEN
div_scal_ser(GEN x,GEN y)2307 div_scal_ser(GEN x, GEN y)
2308 { pari_sp av = avma; return gerepileupto(av, gmul(x, ser_inv(y))); }
2309 static GEN
div_scal_T(GEN x,GEN y,long ty)2310 div_scal_T(GEN x, GEN y, long ty) {
2311 switch(ty)
2312 {
2313 case t_POL: return div_scal_pol(x, y);
2314 case t_SER: return div_scal_ser(x, y);
2315 case t_RFRAC: return div_scal_rfrac(x, y);
2316 }
2317 pari_err_TYPE2("/",x,y);
2318 return NULL; /* LCOV_EXCL_LINE */
2319 }
2320
2321 static GEN
ser2pol_ii(GEN x,long lx,long y1)2322 ser2pol_ii(GEN x, long lx, long y1)
2323 {
2324 GEN y = cgetg(lx, t_POL);
2325 long i;
2326 for (i = lx-1; i > 1; i--) gel(y,i) = gel(x,i);
2327 y[1] = y1; return y;
2328 }
2329
2330 /* assume tx = ty = t_SER, same variable vx */
2331 static GEN
div_ser(GEN x,GEN y,long vx)2332 div_ser(GEN x, GEN y, long vx)
2333 {
2334 long v = valp(x) - valp(y), lx = lg(x), ly = lg(y);
2335 GEN y_lead, z;
2336 pari_sp av = avma;
2337
2338 if (!signe(y)) pari_err_INV("div_ser", y);
2339 if (ser_isexactzero(x))
2340 {
2341 if (lx == 2) return zeroser(vx, v);
2342 return scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), v);
2343 }
2344 y_lead = gel(y,2);
2345 if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
2346 {
2347 GEN y0 = y;
2348 pari_warn(warner,"normalizing a series with 0 leading term");
2349 for (v--, ly--,y++; ly > 2; v--, ly--, y++)
2350 {
2351 y_lead = gel(y,2);
2352 if (!gequal0(y_lead)) break;
2353 }
2354 if (ly <= 2) pari_err_INV("div_ser", y0);
2355 }
2356 if (ly < lx) lx = ly;
2357 z = cgetg(lx,t_SER); z[1] = evalvalp(v) | evalvarn(vx) | evalsigne(1);
2358 x = ser2pol_i(x, lx);
2359 y = ser2pol_ii(y, lx, z[1] & ~VALPBITS);
2360 if (lx == 3) gel(z,2) = gdiv(gel(x,2), gel(y,2));
2361 else
2362 { /* FIXME: use Karp/Markstein */
2363 y = RgXn_mul(x, RgXn_inv_i(y, lx-2), lx-2);
2364 z = fill_ser(z,y);
2365 }
2366 return gerepilecopy(av, z);
2367 }
2368 /* x,y compatible PADIC */
2369 static GEN
divpp(GEN x,GEN y)2370 divpp(GEN x, GEN y) {
2371 pari_sp av;
2372 long a, b;
2373 GEN z, M;
2374
2375 if (!signe(gel(y,4))) pari_err_INV("divpp",y);
2376 if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
2377 a = precp(x);
2378 b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
2379 z = cgetg(5, t_PADIC);
2380 z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
2381 gel(z,2) = icopy(gel(x,2));
2382 gel(z,3) = icopy(M); av = avma;
2383 gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
2384 return z;
2385 }
2386 static GEN
div_polmod_same(GEN T,GEN x,GEN y)2387 div_polmod_same(GEN T, GEN x, GEN y)
2388 {
2389 long v = varn(T);
2390 GEN a, z = cgetg(3, t_POLMOD);
2391 gel(z,1) = RgX_copy(T);
2392 if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
2393 a = gdiv(x, y);
2394 else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
2395 {
2396 pari_sp av = avma;
2397 a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
2398 }
2399 else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
2400 {
2401 pari_sp av = avma;
2402 a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
2403 a = RgX_Rg_div(a, quad_polmod_norm(y, T));
2404 a = gerepileupto(av, a);
2405 }
2406 else
2407 {
2408 pari_sp av = avma;
2409 a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
2410 a = gerepileupto(av, a);
2411 }
2412 gel(z,2) = a; return z;
2413 }
2414 GEN
gdiv(GEN x,GEN y)2415 gdiv(GEN x, GEN y)
2416 {
2417 long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
2418 pari_sp av, tetpil;
2419 GEN z, p1, p2;
2420
2421 if (tx == ty) switch(tx)
2422 {
2423 case t_INT:
2424 return Qdivii(x,y);
2425
2426 case t_REAL: return divrr(x,y);
2427 case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
2428 z = cgetg(3,t_INTMOD);
2429 if (X==Y || equalii(X,Y))
2430 return div_intmod_same(z, X, gel(x,2), gel(y,2));
2431 gel(z,1) = gcdii(X,Y);
2432 warn_coercion(X,Y,gel(z,1));
2433 av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
2434 gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
2435 }
2436 case t_FRAC: {
2437 GEN x1 = gel(x,1), x2 = gel(x,2);
2438 GEN y1 = gel(y,1), y2 = gel(y,2);
2439 z = cgetg(3, t_FRAC);
2440 p1 = gcdii(x1, y1);
2441 if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
2442 p1 = gcdii(x2, y2);
2443 if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
2444 tetpil = avma;
2445 gel(z,2) = mulii(x2,y1);
2446 gel(z,1) = mulii(x1,y2);
2447 normalize_frac(z);
2448 fix_frac_if_int_GC(z,tetpil);
2449 return z;
2450 }
2451 case t_COMPLEX:
2452 if (isintzero(gel(y,1)))
2453 {
2454 y = gel(y,2);
2455 if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
2456 z = cgetg(3,t_COMPLEX);
2457 gel(z,1) = gdiv(gel(x,2), y);
2458 av = avma;
2459 gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
2460 return z;
2461 }
2462 av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
2463 return gerepile(av, tetpil, gdiv(p2,p1));
2464
2465 case t_PADIC:
2466 if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
2467 return divpp(x, y);
2468
2469 case t_QUAD:
2470 if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
2471 av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
2472 return gerepile(av, tetpil, gdiv(p2,p1));
2473
2474 case t_FFELT: return FF_div(x,y);
2475
2476 case t_POLMOD:
2477 if (RgX_equal_var(gel(x,1), gel(y,1)))
2478 z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
2479 else {
2480 av = avma;
2481 z = gerepileupto(av, gmul(x, ginv(y)));
2482 }
2483 return z;
2484
2485 case t_POL:
2486 vx = varn(x);
2487 vy = varn(y);
2488 if (vx != vy) {
2489 if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2490 else return div_scal_pol(x, y);
2491 }
2492 if (!signe(y)) pari_err_INV("gdiv",y);
2493 if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2494 av = avma;
2495 return gerepileupto(av, gred_rfrac2(x,y));
2496
2497 case t_SER:
2498 vx = varn(x);
2499 vy = varn(y);
2500 if (vx != vy) {
2501 if (varncmp(vx, vy) < 0)
2502 {
2503 if (!signe(y)) pari_err_INV("gdiv",y);
2504 return div_ser_scal(x, y);
2505 }
2506 return div_scal_ser(x, y);
2507 }
2508 return div_ser(x, y, vx);
2509 case t_RFRAC:
2510 vx = varn(gel(x,2));
2511 vy = varn(gel(y,2));
2512 if (vx != vy) {
2513 if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2514 else return div_scal_rfrac(x, y);
2515 }
2516 return div_rfrac(x,y);
2517
2518 case t_QFI: av = avma; return gerepileupto(av, qficomp(x, ginv(y)));
2519 case t_QFR: av = avma; return gerepileupto(av, qfrcomp(x, ginv(y)));
2520
2521 case t_MAT:
2522 av = avma; p1 = RgM_inv(y);
2523 if (!p1) pari_err_INV("gdiv",y);
2524 return gerepileupto(av, RgM_mul(x, p1));
2525
2526 default: pari_err_TYPE2("/",x,y);
2527 }
2528
2529 if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2530 {
2531 long s = signe(x);
2532 if (!s) {
2533 if (gequal0(y)) pari_err_INV("gdiv",y);
2534 switch (ty)
2535 {
2536 default: return gen_0;
2537 case t_INTMOD:
2538 z = cgetg(3,t_INTMOD);
2539 gel(z,1) = icopy(gel(y,1));
2540 gel(z,2) = gen_0; return z;
2541 case t_FFELT: return FF_zero(y);
2542 }
2543 }
2544 if (is_pm1(x)) {
2545 if (s > 0) return ginv(y);
2546 av = avma; return gerepileupto(av, ginv(gneg(y)));
2547 }
2548 switch(ty)
2549 {
2550 case t_REAL: return divir(x,y);
2551 case t_INTMOD:
2552 z = cgetg(3, t_INTMOD);
2553 return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2554 case t_FRAC:
2555 z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2556 if (equali1(p1))
2557 {
2558 set_avma((pari_sp)z);
2559 gel(z,2) = icopy(gel(y,1));
2560 gel(z,1) = mulii(gel(y,2), x);
2561 normalize_frac(z);
2562 fix_frac_if_int(z);
2563 }
2564 else
2565 {
2566 x = diviiexact(x,p1); tetpil = avma;
2567 gel(z,2) = diviiexact(gel(y,1), p1);
2568 gel(z,1) = mulii(gel(y,2), x);
2569 normalize_frac(z);
2570 fix_frac_if_int_GC(z,tetpil);
2571 }
2572 return z;
2573
2574 case t_FFELT: return Z_FF_div(x,y);
2575 case t_COMPLEX: return divRc(x,y);
2576 case t_PADIC: return divTp(x, y);
2577 case t_QUAD:
2578 av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
2579 return gerepile(av, tetpil, gdiv(p2,p1));
2580 }
2581 }
2582 if (gequal0(y))
2583 {
2584 if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2585 if (ty != t_MAT) pari_err_INV("gdiv",y);
2586 }
2587
2588 if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2589 {
2590 case t_REAL:
2591 switch(ty)
2592 {
2593 case t_INT: return divri(x,y);
2594 case t_FRAC:
2595 av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2596 return gerepileuptoleaf(av, z);
2597 case t_COMPLEX: return divRc(x, y);
2598 case t_QUAD: return divfq(x, y, realprec(x));
2599 default: pari_err_TYPE2("/",x,y);
2600 }
2601
2602 case t_INTMOD:
2603 switch(ty)
2604 {
2605 case t_INT:
2606 z = cgetg(3, t_INTMOD);
2607 return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2608 case t_FRAC: { GEN X = gel(x,1);
2609 z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2610 return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2611 }
2612 case t_FFELT:
2613 if (!equalii(gel(x,1),FF_p_i(y)))
2614 pari_err_OP("/",x,y);
2615 return Z_FF_div(gel(x,2),y);
2616
2617 case t_COMPLEX:
2618 av = avma;
2619 return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
2620
2621 case t_QUAD:
2622 av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
2623 return gerepile(av,tetpil, gdiv(p2,p1));
2624
2625 case t_PADIC: { GEN X = gel(x,1);
2626 z = cgetg(3, t_INTMOD);
2627 return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2628 }
2629 case t_REAL: pari_err_TYPE2("/",x,y);
2630 }
2631
2632 case t_FRAC:
2633 switch(ty)
2634 {
2635 case t_INT: z = cgetg(3, t_FRAC);
2636 p1 = gcdii(y,gel(x,1));
2637 if (equali1(p1))
2638 {
2639 set_avma((pari_sp)z); tetpil = 0;
2640 gel(z,1) = icopy(gel(x,1));
2641 }
2642 else
2643 {
2644 y = diviiexact(y,p1); tetpil = avma;
2645 gel(z,1) = diviiexact(gel(x,1), p1);
2646 }
2647 gel(z,2) = mulii(gel(x,2),y);
2648 normalize_frac(z);
2649 if (tetpil) fix_frac_if_int_GC(z,tetpil);
2650 return z;
2651
2652 case t_REAL:
2653 av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
2654 return gerepile(av, tetpil, divir(gel(x,1), p1));
2655
2656 case t_INTMOD: { GEN Y = gel(y,1);
2657 z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2658 return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
2659 }
2660
2661 case t_FFELT: av=avma;
2662 return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
2663
2664 case t_COMPLEX: return divRc(x, y);
2665
2666 case t_PADIC:
2667 if (!signe(gel(x,1))) return gen_0;
2668 return divTp(x, y);
2669
2670 case t_QUAD:
2671 av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
2672 return gerepile(av,tetpil,gdiv(p2,p1));
2673 }
2674
2675 case t_FFELT:
2676 switch (ty)
2677 {
2678 case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
2679 case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
2680 case t_INTMOD:
2681 if (!equalii(gel(y,1),FF_p_i(x)))
2682 pari_err_OP("/",x,y);
2683 return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
2684 default:
2685 pari_err_TYPE2("/",x,y);
2686 }
2687 break;
2688
2689 case t_COMPLEX:
2690 switch(ty)
2691 {
2692 case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
2693 case t_INTMOD: return mulRc_direct(ginv(y), x);
2694 case t_PADIC:
2695 return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
2696 case t_QUAD:
2697 lx = precision(x); if (!lx) pari_err_OP("/",x,y);
2698 return divfq(x, y, lx);
2699 }
2700
2701 case t_PADIC:
2702 switch(ty)
2703 {
2704 case t_INT: case t_FRAC: { GEN p = gel(x,2);
2705 return signe(gel(x,4))? divpT(x, y)
2706 : zeropadic(p, valp(x) - Q_pval(y,p));
2707 }
2708 case t_INTMOD: { GEN Y = gel(y,1);
2709 z = cgetg(3, t_INTMOD);
2710 return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2711 }
2712 case t_COMPLEX: case t_QUAD:
2713 av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
2714 return gerepile(av,tetpil,gdiv(p1,p2));
2715
2716 case t_REAL: pari_err_TYPE2("/",x,y);
2717 }
2718
2719 case t_QUAD:
2720 switch (ty)
2721 {
2722 case t_INT: case t_INTMOD: case t_FRAC:
2723 z = cgetg(4,t_QUAD);
2724 gel(z,1) = ZX_copy(gel(x,1));
2725 gel(z,2) = gdiv(gel(x,2), y);
2726 gel(z,3) = gdiv(gel(x,3), y); return z;
2727 case t_REAL: return divqf(x, y, realprec(y));
2728 case t_PADIC: return divTp(x, y);
2729 case t_COMPLEX:
2730 ly = precision(y); if (!ly) pari_err_OP("/",x,y);
2731 return divqf(x, y, ly);
2732 }
2733 }
2734 switch(ty) {
2735 case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2736 return gmul(x, ginv(y)); /* missing gerepile, for speed */
2737 case t_MAT:
2738 av = avma; p1 = RgM_inv(y);
2739 if (!p1) pari_err_INV("gdiv",y);
2740 return gerepileupto(av, gmul(x, p1));
2741 case t_VEC: case t_COL:
2742 case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2743 pari_err_TYPE2("/",x,y);
2744 }
2745 switch(tx) {
2746 case t_VEC: case t_COL: case t_MAT:
2747 z = cgetg_copy(x, &lx);
2748 for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
2749 return z;
2750 case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
2751 pari_err_TYPE2("/",x,y);
2752 }
2753
2754 vy = gvar(y);
2755 if (tx == t_POLMOD) { GEN X = gel(x,1);
2756 vx = varn(X);
2757 if (vx != vy) {
2758 if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2759 retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
2760 }
2761 /* y is POL, SER or RFRAC */
2762 av = avma;
2763 switch(ty)
2764 {
2765 case t_RFRAC: y = gmod(ginv(y), X); break;
2766 default: y = ginvmod(gmod(y,X), X);
2767 }
2768 return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
2769 }
2770 /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2771 * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
2772 * variable NO_VARIABLE, then the operation is incorrect */
2773 vx = gvar(x);
2774 if (vx != vy) { /* includes cases where one is scalar */
2775 if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2776 else return div_scal_T(x, y, ty);
2777 }
2778 switch(tx)
2779 {
2780 case t_POL:
2781 switch(ty)
2782 {
2783 case t_SER:
2784 if (lg(y) == 2)
2785 return zeroser(vx, RgX_val(x) - valp(y));
2786 p1 = RgX_to_ser(x,lg(y));
2787 p2 = div_ser(p1, y, vx);
2788 settyp(p1, t_VECSMALL); /* p1 left on stack */
2789 return p2;
2790
2791 case t_RFRAC:
2792 {
2793 GEN y1 = gel(y,1), y2 = gel(y,2);
2794 if (typ(y1) == t_POL && varn(y1) == vx)
2795 return mul_rfrac_scal(y2, y1, x);
2796 av = avma;
2797 return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2798 }
2799 }
2800 break;
2801
2802 case t_SER:
2803 switch(ty)
2804 {
2805 case t_POL:
2806 if (lg(x) == 2)
2807 return zeroser(vx, valp(x) - RgX_val(y));
2808 p1 = RgX_to_ser_inexact(y,lg(x));
2809 p2 = div_ser(x, p1, vx);
2810 settyp(p1, t_VECSMALL); /* p1 left on stack */
2811 return p2;
2812 case t_RFRAC:
2813 av = avma;
2814 return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2815 }
2816 break;
2817
2818 case t_RFRAC:
2819 switch(ty)
2820 {
2821 case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2822 case t_SER:
2823 av = avma; z = RgX_to_ser_inexact(gel(x,2), lg(y));
2824 return gerepileupto(av, gdiv(gel(x,1), gmul(z,y)));
2825 }
2826 break;
2827 }
2828 pari_err_TYPE2("/",x,y);
2829 return NULL; /* LCOV_EXCL_LINE */
2830 }
2831
2832 /********************************************************************/
2833 /** **/
2834 /** SIMPLE MULTIPLICATION **/
2835 /** **/
2836 /********************************************************************/
2837 GEN
gmulsg(long s,GEN y)2838 gmulsg(long s, GEN y)
2839 {
2840 long ly, i;
2841 pari_sp av;
2842 GEN z;
2843
2844 switch(typ(y))
2845 {
2846 case t_INT: return mulsi(s,y);
2847 case t_REAL: return s? mulsr(s,y): gen_0; /* gmul semantic */
2848 case t_INTMOD: { GEN p = gel(y,1);
2849 z = cgetg(3,t_INTMOD);
2850 gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
2851 gel(z,1) = icopy(p); return z;
2852 }
2853 case t_FFELT: return FF_Z_mul(y,stoi(s));
2854 case t_FRAC:
2855 if (!s) return gen_0;
2856 z = cgetg(3,t_FRAC);
2857 i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
2858 if (i == 1)
2859 {
2860 gel(z,2) = icopy(gel(y,2));
2861 gel(z,1) = mulis(gel(y,1), s);
2862 }
2863 else
2864 {
2865 gel(z,2) = divis(gel(y,2), i);
2866 gel(z,1) = mulis(gel(y,1), s/i);
2867 fix_frac_if_int(z);
2868 }
2869 return z;
2870
2871 case t_COMPLEX:
2872 if (!s) return gen_0;
2873 z = cgetg(3, t_COMPLEX);
2874 gel(z,1) = gmulsg(s,gel(y,1));
2875 gel(z,2) = gmulsg(s,gel(y,2)); return z;
2876
2877 case t_PADIC:
2878 if (!s) return gen_0;
2879 av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
2880
2881 case t_QUAD: z = cgetg(4, t_QUAD);
2882 gel(z,1) = ZX_copy(gel(y,1));
2883 gel(z,2) = gmulsg(s,gel(y,2));
2884 gel(z,3) = gmulsg(s,gel(y,3)); return z;
2885
2886 case t_POLMOD:
2887 retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
2888
2889 case t_POL:
2890 if (!signe(y)) return RgX_copy(y);
2891 if (!s) return scalarpol(Rg_get_0(y), varn(y));
2892 z = cgetg_copy(y, &ly); z[1]=y[1];
2893 for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2894 return normalizepol_lg(z, ly);
2895
2896 case t_SER:
2897 if (ser_isexactzero(y)) return gcopy(y);
2898 if (!s) return Rg_get_0(y);
2899 z = cgetg_copy(y, &ly); z[1]=y[1];
2900 for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2901 return normalize(z);
2902
2903 case t_RFRAC:
2904 if (!s) return zeropol(varn(gel(y,2)));
2905 if (s == 1) return gcopy(y);
2906 if (s == -1) return gneg(y);
2907 return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
2908
2909 case t_VEC: case t_COL: case t_MAT:
2910 z = cgetg_copy(y, &ly);
2911 for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2912 return z;
2913 }
2914 pari_err_TYPE("gmulsg",y);
2915 return NULL; /* LCOV_EXCL_LINE */
2916 }
2917
2918 /********************************************************************/
2919 /** **/
2920 /** SIMPLE DIVISION **/
2921 /** **/
2922 /********************************************************************/
2923
2924 GEN
gdivgs(GEN x,long s)2925 gdivgs(GEN x, long s)
2926 {
2927 long tx = typ(x), lx, i;
2928 pari_sp av;
2929 GEN z;
2930
2931 if (!s)
2932 {
2933 if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
2934 pari_err_INV("gdivgs",gen_0);
2935 }
2936 switch(tx)
2937 {
2938 case t_INT: return Qdivis(x, s);
2939 case t_REAL: return divrs(x,s);
2940
2941 case t_INTMOD:
2942 z = cgetg(3, t_INTMOD);
2943 return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
2944
2945 case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
2946
2947 case t_FRAC: z = cgetg(3, t_FRAC);
2948 i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
2949 if (i == 1)
2950 {
2951 gel(z,2) = mulsi(s, gel(x,2));
2952 gel(z,1) = icopy(gel(x,1));
2953 }
2954 else
2955 {
2956 gel(z,2) = mulsi(s/i, gel(x,2));
2957 gel(z,1) = divis(gel(x,1), i);
2958 }
2959 normalize_frac(z);
2960 fix_frac_if_int(z); return z;
2961
2962 case t_COMPLEX: z = cgetg(3, t_COMPLEX);
2963 gel(z,1) = gdivgs(gel(x,1),s);
2964 gel(z,2) = gdivgs(gel(x,2),s); return z;
2965
2966 case t_PADIC: /* divpT */
2967 {
2968 GEN p = gel(x,2);
2969 if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
2970 av = avma;
2971 return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
2972 }
2973
2974 case t_QUAD: z = cgetg(4, t_QUAD);
2975 gel(z,1) = ZX_copy(gel(x,1));
2976 gel(z,2) = gdivgs(gel(x,2),s);
2977 gel(z,3) = gdivgs(gel(x,3),s); return z;
2978
2979 case t_POLMOD:
2980 retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
2981
2982 case t_RFRAC:
2983 if (s == 1) return gcopy(x);
2984 else if (s == -1) return gneg(x);
2985 return div_rfrac_scal(x, stoi(s));
2986
2987 case t_POL: case t_SER:
2988 z = cgetg_copy(x, &lx); z[1] = x[1];
2989 for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
2990 return z;
2991 case t_VEC: case t_COL: case t_MAT:
2992 z = cgetg_copy(x, &lx);
2993 for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
2994 return z;
2995
2996 }
2997 pari_err_TYPE2("/",x, stoi(s));
2998 return NULL; /* LCOV_EXCL_LINE */
2999 }
3000
3001 /* True shift (exact multiplication by 2^n) */
3002 GEN
gmul2n(GEN x,long n)3003 gmul2n(GEN x, long n)
3004 {
3005 long lx, i, k, l;
3006 GEN z, a, b;
3007
3008 switch(typ(x))
3009 {
3010 case t_INT:
3011 if (n>=0) return shifti(x,n);
3012 if (!signe(x)) return gen_0;
3013 l = vali(x); n = -n;
3014 if (n<=l) return shifti(x,-n);
3015 z = cgetg(3,t_FRAC);
3016 gel(z,1) = shifti(x,-l);
3017 gel(z,2) = int2n(n-l); return z;
3018
3019 case t_REAL:
3020 return shiftr(x,n);
3021
3022 case t_INTMOD: b = gel(x,1); a = gel(x,2);
3023 z = cgetg(3,t_INTMOD);
3024 if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
3025 gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
3026 gel(z,1) = icopy(b); return z;
3027
3028 case t_FFELT: return FF_mul2n(x,n);
3029
3030 case t_FRAC: a = gel(x,1); b = gel(x,2);
3031 l = vali(a);
3032 k = vali(b);
3033 if (n+l >= k)
3034 {
3035 if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
3036 l = n-k; k = -k;
3037 }
3038 else
3039 {
3040 k = -(l+n); l = -l;
3041 }
3042 z = cgetg(3,t_FRAC);
3043 gel(z,1) = shifti(a,l);
3044 gel(z,2) = shifti(b,k); return z;
3045
3046 case t_COMPLEX: z = cgetg(3,t_COMPLEX);
3047 gel(z,1) = gmul2n(gel(x,1),n);
3048 gel(z,2) = gmul2n(gel(x,2),n); return z;
3049
3050 case t_QUAD: z = cgetg(4,t_QUAD);
3051 gel(z,1) = ZX_copy(gel(x,1));
3052 gel(z,2) = gmul2n(gel(x,2),n);
3053 gel(z,3) = gmul2n(gel(x,3),n); return z;
3054
3055 case t_POLMOD:
3056 retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
3057
3058 case t_POL:
3059 z = cgetg_copy(x, &lx); z[1] = x[1];
3060 for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3061 return normalizepol_lg(z, lx); /* needed if char = 2 */
3062 case t_SER:
3063 if (ser_isexactzero(x)) return gcopy(x);
3064 z = cgetg_copy(x, &lx); z[1] = x[1];
3065 for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3066 return normalize(z); /* needed if char = 2 */
3067 case t_VEC: case t_COL: case t_MAT:
3068 z = cgetg_copy(x, &lx);
3069 for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
3070 return z;
3071
3072 case t_RFRAC: /* int2n wrong if n < 0 */
3073 return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
3074
3075 case t_PADIC: /* int2n wrong if n < 0 */
3076 return gmul(gmul2n(gen_1,n),x);
3077 }
3078 pari_err_TYPE("gmul2n",x);
3079 return NULL; /* LCOV_EXCL_LINE */
3080 }
3081
3082 /*******************************************************************/
3083 /* */
3084 /* INVERSE */
3085 /* */
3086 /*******************************************************************/
3087 static GEN
inv_polmod(GEN T,GEN x)3088 inv_polmod(GEN T, GEN x)
3089 {
3090 GEN z = cgetg(3,t_POLMOD), a;
3091 gel(z,1) = RgX_copy(T);
3092 if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
3093 a = ginv(x);
3094 else
3095 {
3096 if (lg(T) == 5) /* quadratic fields */
3097 a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
3098 else
3099 a = RgXQ_inv(x, T);
3100 }
3101 gel(z,2) = a; return z;
3102 }
3103
3104 GEN
ginv(GEN x)3105 ginv(GEN x)
3106 {
3107 long s;
3108 pari_sp av, tetpil;
3109 GEN z, y, p1, p2;
3110
3111 switch(typ(x))
3112 {
3113 case t_INT:
3114 if (is_pm1(x)) return icopy(x);
3115 s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
3116 z = cgetg(3,t_FRAC);
3117 gel(z,1) = s<0? gen_m1: gen_1;
3118 gel(z,2) = absi(x); return z;
3119
3120 case t_REAL: return invr(x);
3121
3122 case t_INTMOD: z=cgetg(3,t_INTMOD);
3123 gel(z,1) = icopy(gel(x,1));
3124 gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
3125
3126 case t_FRAC: {
3127 GEN a = gel(x,1), b = gel(x,2);
3128 s = signe(a);
3129 if (is_pm1(a)) return s > 0? icopy(b): negi(b);
3130 z = cgetg(3,t_FRAC);
3131 gel(z,1) = icopy(b);
3132 gel(z,2) = icopy(a);
3133 normalize_frac(z); return z;
3134 }
3135 case t_COMPLEX:
3136 av=avma;
3137 p1=cxnorm(x);
3138 p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
3139 tetpil=avma;
3140 return gerepile(av,tetpil,divcR(p2,p1));
3141
3142 case t_QUAD:
3143 av=avma; p1=gnorm(x); p2=conj_i(x); tetpil=avma;
3144 return gerepile(av,tetpil,gdiv(p2,p1));
3145
3146 case t_PADIC: z = cgetg(5,t_PADIC);
3147 if (!signe(gel(x,4))) pari_err_INV("ginv",x);
3148 z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
3149 gel(z,2) = icopy(gel(x,2));
3150 gel(z,3) = icopy(gel(x,3));
3151 gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
3152
3153 case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
3154 case t_FFELT: return FF_inv(x);
3155 case t_POL: return gred_rfrac_simple(gen_1,x);
3156 case t_SER: return ser_inv(x);
3157 case t_RFRAC:
3158 {
3159 GEN n = gel(x,1), d = gel(x,2);
3160 pari_sp av = avma, ltop;
3161 if (gequal0(n)) pari_err_INV("ginv",x);
3162
3163 n = simplify_shallow(n);
3164 if (typ(n) != t_POL || varn(n) != varn(d))
3165 {
3166 if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
3167 ltop = avma;
3168 z = RgX_Rg_div(d,n);
3169 } else {
3170 ltop = avma;
3171 z = cgetg(3,t_RFRAC);
3172 gel(z,1) = RgX_copy(d);
3173 gel(z,2) = RgX_copy(n);
3174 }
3175 stackdummy(av, ltop);
3176 return z;
3177 }
3178
3179 case t_QFR:
3180 av = avma; z = cgetg(5, t_QFR);
3181 gel(z,1) = gel(x,1);
3182 gel(z,2) = negi( gel(x,2) );
3183 gel(z,3) = gel(x,3);
3184 gel(z,4) = negr( gel(x,4) );
3185 return gerepileupto(av, redreal(z));
3186
3187 case t_QFI:
3188 y = gcopy(x);
3189 if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
3190 togglesign(gel(y,2));
3191 return y;
3192 case t_MAT:
3193 y = RgM_inv(x);
3194 if (!y) pari_err_INV("ginv",x);
3195 return y;
3196 case t_VECSMALL:
3197 {
3198 long i, lx = lg(x)-1;
3199 y = zero_zv(lx);
3200 for (i=1; i<=lx; i++)
3201 {
3202 long xi = x[i];
3203 if (xi<1 || xi>lx || y[xi])
3204 pari_err_TYPE("ginv [not a permutation]", x);
3205 y[xi] = i;
3206 }
3207 return y;
3208 }
3209 }
3210 pari_err_TYPE("inverse",x);
3211 return NULL; /* LCOV_EXCL_LINE */
3212 }
3213