1 /* $Id: gen1.c 11470 2008-12-19 19:33:27Z kb $
2
3 Copyright (C) 2000 The PARI group.
4
5 This file is part of the PARI/GP package.
6
7 PARI/GP is free software; you can redistribute it and/or modify it under the
8 terms of the GNU General Public License as published by the Free Software
9 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15
16 /********************************************************************/
17 /** **/
18 /** GENERIC OPERATIONS **/
19 /** (first part) **/
20 /** **/
21 /********************************************************************/
22 #include "pari.h"
23 #include "paripriv.h"
24
25 #define fix_frac(z) if (signe(z[2])<0)\
26 {\
27 setsigne(z[1],-signe(z[1]));\
28 setsigne(z[2],1);\
29 }
30
31 /* assume z[1] was created last */
32 #define fix_frac_if_int(z) if (is_pm1(z[2]))\
33 z = gerepileupto((pari_sp)(z+3), gel(z,1));
34
35 /* assume z[1] was created last */
36 #define fix_frac_if_int_GC(z,tetpil) { if (is_pm1(z[2]))\
37 z = gerepileupto((pari_sp)(z+3), gel(z,1));\
38 else\
39 gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
40
41 static long
kro_quad(GEN x,GEN y)42 kro_quad(GEN x, GEN y)
43 {
44 long k;
45 pari_sp av=avma;
46
47 x = subii(sqri(gel(x,3)), shifti(gel(x,2),2));
48 k = kronecker(x,y); avma=av; return k;
49 }
50
51 /* is -1 not a square in Zp, assume p prime */
52 INLINE int
Zp_nosquare_m1(GEN p)53 Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
54
55 static GEN addpp(GEN x, GEN y);
56 static GEN mulpp(GEN x, GEN y);
57 static GEN divpp(GEN x, GEN y);
58 /* Argument codes for inline routines
59 * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
60 * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
61 * T: some type (to be converted to PADIC)
62 */
63 static GEN
addRc(GEN x,GEN y)64 addRc(GEN x, GEN y) {
65 GEN z = cgetg(3,t_COMPLEX);
66 gel(z,1) = gadd(x,gel(y,1));
67 gel(z,2) = gcopy(gel(y,2)); return z;
68 }
69 static GEN
mulRc(GEN x,GEN y)70 mulRc(GEN x, GEN y) {
71 GEN z = cgetg(3,t_COMPLEX);
72 gel(z,1) = gmul(x,gel(y,1));
73 gel(z,2) = gmul(x,gel(y,2)); return z;
74 }
75 static GEN
divRc(GEN x,GEN y)76 divRc(GEN x, GEN y) {
77 GEN a, b, N, z = cgetg(3,t_COMPLEX);
78 pari_sp tetpil, av = avma;
79 a = gmul(x, gel(y,1));
80 b = gmul(x, gel(y,2)); if(!gcmp0(b)) b = gneg_i(b);
81 N = cxnorm(y); tetpil = avma;
82 gel(z,1) = gdiv(a, N);
83 gel(z,2) = gdiv(b, N); gerepilecoeffssp(av,tetpil,z+1,2); return z;
84 }
85 static GEN
divcR(GEN x,GEN y)86 divcR(GEN x, GEN y) {
87 GEN z = cgetg(3,t_COMPLEX);
88 gel(z,1) = gdiv(gel(x,1), y);
89 gel(z,2) = gdiv(gel(x,2), y); return z;
90 }
91 static GEN
addRq(GEN x,GEN y)92 addRq(GEN x, GEN y) {
93 GEN z = cgetg(4,t_QUAD);
94 gel(z,1) = gcopy(gel(y,1));
95 gel(z,2) = gadd(x, gel(y,2));
96 gel(z,3) = gcopy(gel(y,3)); return z;
97 }
98 static GEN
mulRq(GEN x,GEN y)99 mulRq(GEN x, GEN y) {
100 GEN z = cgetg(4,t_QUAD);
101 gel(z,1) = gcopy(gel(y,1));
102 gel(z,2) = gmul(x,gel(y,2));
103 gel(z,3) = gmul(x,gel(y,3)); return z;
104 }
105 static GEN
addqf(GEN x,GEN y,long prec)106 addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
107 long i = gexpo(x) - gexpo(y);
108 if (i <= 0) i = 0; else i >>= TWOPOTBITS_IN_LONG;
109 return gerepileupto(av, gadd(y, quadtoc(x, prec + i)));
110 }
111 static GEN
mulqf(GEN x,GEN y,long prec)112 mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
113 return gerepileupto(av, gmul(y, quadtoc(x, prec)));
114 }
115 static GEN
divqf(GEN x,GEN y,long prec)116 divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
117 return gerepileupto(av, gdiv(quadtoc(x,prec), y));
118 }
119 static GEN
divfq(GEN x,GEN y,long prec)120 divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
121 return gerepileupto(av, gdiv(x, quadtoc(y,prec)));
122 }
123 /* y PADIC, x + y by converting x to padic */
124 static GEN
addTp(GEN x,GEN y)125 addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
126
127 if (!valp(y)) z = cvtop2(x,y);
128 else {
129 long l = signe(y[4])? valp(y) + precp(y): valp(y);
130 z = cvtop(x, gel(y,2), l);
131 }
132 return gerepileupto(av, addpp(z, y));
133 }
134 /* y PADIC, x * y by converting x to padic */
135 static GEN
mulTp(GEN x,GEN y)136 mulTp(GEN x, GEN y) { pari_sp av = avma;
137 return gerepileupto(av, mulpp(cvtop2(x,y), y));
138 }
139 /* y PADIC, non zero x / y by converting x to padic */
140 static GEN
divTp(GEN x,GEN y)141 divTp(GEN x, GEN y) { pari_sp av = avma;
142 return gerepileupto(av, divpp(cvtop2(x,y), y));
143 }
144 /* x PADIC, x / y by converting y to padic */
145 static GEN
divpT(GEN x,GEN y)146 divpT(GEN x, GEN y) { pari_sp av = avma;
147 return gerepileupto(av, divpp(x, cvtop2(y,x)));
148 }
149
150 /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
151 * clean memory from z on */
152 static GEN
add_intmod_same(GEN z,GEN X,GEN x,GEN y)153 add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
154 if (lgefint(X) == 3) {
155 ulong u = Fl_add(itou(x),itou(y), X[2]);
156 avma = (pari_sp)z; gel(z,2) = utoi(u);
157 }
158 else {
159 GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
160 gel(z,2) = gerepileuptoint((pari_sp)z, u);
161 }
162 gel(z,1) = icopy(X); return z;
163 }
164 /* cf add_intmod_same */
165 static GEN
mul_intmod_same(GEN z,GEN X,GEN x,GEN y)166 mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
167 if (lgefint(X) == 3) {
168 ulong u = Fl_mul(itou(x),itou(y), X[2]);
169 avma = (pari_sp)z; gel(z,2) = utoi(u);
170 }
171 else
172 gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
173 gel(z,1) = icopy(X); return z;
174 }
175 /* cf add_intmod_same */
176 static GEN
div_intmod_same(GEN z,GEN X,GEN x,GEN y)177 div_intmod_same(GEN z, GEN X, GEN x, GEN y)
178 {
179 if (lgefint(X) == 3) {
180 ulong m = (ulong)X[2], u = Fl_div(itou(x), itou(y), m);
181 avma = (pari_sp)z; gel(z,2) = utoi(u);
182 }
183 else
184 gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
185 gel(z,1) = icopy(X); return z;
186 }
187
188 /*******************************************************************/
189 /* */
190 /* REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC) */
191 /* */
192 /* (static routines are not memory clean, but OK for gerepileupto) */
193 /*******************************************************************/
194 /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
195 GEN
gred_rfrac_simple(GEN n,GEN d)196 gred_rfrac_simple(GEN n, GEN d)
197 {
198 GEN c, cn, cd, z;
199
200 cd = content(d);
201 cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
202 if (!gcmp1(cd)) {
203 d = RgX_Rg_div(d,cd);
204 if (!gcmp1(cn))
205 {
206 if (gcmp0(cn)) {
207 n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
208 c = gen_1;
209 } else {
210 n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
211 c = gdiv(cn,cd);
212 }
213 }
214 else
215 c = ginv(cd);
216 } else {
217 if (!gcmp1(cn))
218 {
219 if (gcmp0(cn)) {
220 c = gen_1;
221 } else {
222 n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
223 c = cn;
224 }
225 } else {
226 GEN y = cgetg(3,t_RFRAC);
227 gel(y,1) = gcopy(n);
228 gel(y,2) = gcopy(d); return y;
229 }
230 }
231
232 if (typ(c) == t_POL)
233 {
234 z = c;
235 do { z = content(z); } while (typ(z) == t_POL);
236 cd = denom(z);
237 cn = gmul(c, cd);
238 }
239 else
240 {
241 cn = numer(c);
242 cd = denom(c);
243 }
244 z = cgetg(3,t_RFRAC);
245 gel(z,1) = gmul(n, cn);
246 gel(z,2) = gmul(d, cd); return z;
247 }
248
249 static GEN
fix_rfrac(GEN x,long d)250 fix_rfrac(GEN x, long d)
251 {
252 GEN z, N, D;
253 if (!d) return x;
254 z = cgetg(3, t_RFRAC);
255 N = gel(x,1);
256 D = gel(x,2);
257 if (d > 0) {
258 gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
259 : monomialcopy(N,d,varn(D));
260 gel(z, 2) = gcopy(D);
261 } else {
262 gel(z, 1) = gcopy(N);
263 gel(z, 2) = RgX_shift(D, -d);
264 }
265 return z;
266 }
267
268 /* assume d != 0 */
269 static GEN
gred_rfrac2_i(GEN n,GEN d)270 gred_rfrac2_i(GEN n, GEN d)
271 {
272 GEN y, z;
273 long tn, td, v;
274
275 n = simplify_i(n);
276 if (isexactzero(n)) return gcopy(n);
277 d = simplify_i(d);
278 td = typ(d);
279 if (td!=t_POL || varncmp(varn(d), gvar(n)) > 0) return gdiv(n,d);
280 tn = typ(n);
281 if (tn!=t_POL)
282 {
283 if (varncmp(varn(d), gvar2(n)) < 0) return gred_rfrac_simple(n,d);
284 pari_err(talker,"incompatible variables in gred");
285 }
286 if (varncmp(varn(d), varn(n)) < 0) return gred_rfrac_simple(n,d);
287 if (varncmp(varn(d), varn(n)) > 0) return RgX_Rg_div(n,d);
288
289 /* now n and d are t_POLs in the same variable */
290 v = polvaluation(n, &n) - polvaluation(d, &d);
291 if (!degpol(d))
292 {
293 n = RgX_Rg_div(n,gel(d,2));
294 return v? RgX_mulXn(n,v): n;
295 }
296
297 /* X does not divide gcd(n,d), deg(d) > 0 */
298 if (!isinexact(n) && !isinexact(d))
299 {
300 y = RgX_divrem(n, d, &z);
301 if (!signe(z)) return v? RgX_mulXn(y, v): y;
302 z = srgcd(d, z);
303 if (degpol(z)) { n = gdeuc(n,z); d = gdeuc(d,z); }
304 }
305 return fix_rfrac(gred_rfrac_simple(n,d), v);
306 }
307
308 GEN
gred_rfrac2(GEN x1,GEN x2)309 gred_rfrac2(GEN x1, GEN x2)
310 {
311 pari_sp av = avma;
312 return gerepileupto(av, gred_rfrac2_i(x1, x2));
313 }
314
315 /* x1,x2 t_INT, return x1/x2 in reduced form */
316 GEN
gred_frac2(GEN x1,GEN x2)317 gred_frac2(GEN x1, GEN x2)
318 {
319 GEN p1, y = dvmdii(x1,x2,&p1);
320 pari_sp av;
321
322 if (p1 == gen_0) return y; /* gen_0 intended */
323 av = avma;
324 p1 = gcdii(x2,p1);
325 if (is_pm1(p1))
326 {
327 avma = av; y = cgetg(3,t_FRAC);
328 gel(y,1) = icopy(x1);
329 gel(y,2) = icopy(x2);
330 }
331 else
332 {
333 p1 = gclone(p1);
334 avma = av; y = cgetg(3,t_FRAC);
335 gel(y,1) = diviiexact(x1,p1);
336 gel(y,2) = diviiexact(x2,p1);
337 gunclone(p1);
338 }
339 fix_frac(y); return y;
340 }
341
342 /********************************************************************/
343 /** **/
344 /** SUBTRACTION **/
345 /** **/
346 /********************************************************************/
347
348 GEN
gsub(GEN x,GEN y)349 gsub(GEN x, GEN y)
350 {
351 pari_sp tetpil, av = avma;
352 y = gneg_i(y); tetpil = avma;
353 return gerepile(av,tetpil, gadd(x,y));
354 }
355
356 /********************************************************************/
357 /** **/
358 /** ADDITION **/
359 /** **/
360 /********************************************************************/
361 /* x, y compatible PADIC */
362 static GEN
addpp(GEN x,GEN y)363 addpp(GEN x, GEN y)
364 {
365 pari_sp av = avma;
366 long c,d,e,r,rx,ry;
367 GEN u, z, mod, p = gel(x,2);
368
369 (void)new_chunk(5 + lgefint(x[3]) + lgefint(y[3]));
370 e = valp(x);
371 r = valp(y); d = r-e;
372 if (d < 0) { swap(x,y); e = r; d = -d; }
373 rx = precp(x);
374 ry = precp(y);
375 if (d) /* v(x) < v(y) */
376 {
377 r = d+ry; z = powiu(p,d);
378 if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
379 u = addii(gel(x,4), mulii(z,gel(y,4)));
380 u = remii(u, mod);
381 }
382 else
383 {
384 if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
385 u = addii(gel(x,4), gel(y,4));
386 if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
387 {
388 avma = av; return zeropadic(p, e+r);
389 }
390 if (c)
391 {
392 mod = diviiexact(mod, powiu(p,c));
393 r -= c;
394 e += c;
395 }
396 u = remii(u, mod);
397 }
398 avma = av; z = cgetg(5,t_PADIC);
399 z[1] = evalprecp(r) | evalvalp(e);
400 gel(z,2) = icopy(p);
401 gel(z,3) = icopy(mod);
402 gel(z,4) = icopy(u); return z;
403 }
404
405 /* return x + y, where x is t_INT or t_FRAC(N), y t_PADIC */
406 static GEN
addQp(GEN x,GEN y)407 addQp(GEN x, GEN y)
408 {
409 pari_sp av;
410 long tx,vy,py,d,r,e;
411 GEN z,q,p,p1,p2,mod,u;
412
413 if (gcmp0(x)) return gcopy(y);
414
415 av = avma; p = gel(y,2); tx = typ(x);
416 e = (tx == t_INT)? Z_pvalrem(x,p,&p1)
417 : Z_pvalrem(gel(x,1),p,&p1) -
418 Z_pvalrem(gel(x,2),p,&p2);
419 vy = valp(y); d = vy - e; py = precp(y); r = d + py;
420 if (r <= 0) { avma = av; return gcopy(y); }
421 mod = gel(y,3);
422 u = gel(y,4);
423 (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
424
425 if (d > 0)
426 {
427 q = powiu(p,d);
428 mod = mulii(mod, q);
429 u = mulii(u, q);
430 if (tx != t_INT && !is_pm1(p2)) p1 = mulii(p1, Fp_inv(p2,mod));
431 u = addii(u, p1);
432 }
433 else if (d < 0)
434 {
435 q = powiu(p,-d);
436 if (tx != t_INT && !is_pm1(p2)) p1 = mulii(p1, Fp_inv(p2,mod));
437 p1 = mulii(p1, q);
438 u = addii(u, p1);
439 r = py; e = vy;
440 }
441 else
442 {
443 long c;
444 if (tx != t_INT && !is_pm1(p2)) p1 = mulii(p1, Fp_inv(p2,mod));
445 u = addii(u, p1);
446 if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
447 {
448 avma = av; return zeropadic(p,e+r);
449 }
450 if (c)
451 {
452 mod = diviiexact(mod, powiu(p,c));
453 r -= c;
454 e += c;
455 }
456 }
457 u = modii(u, mod);
458 avma = av; z = cgetg(5,t_PADIC);
459 z[1] = evalprecp(r) | evalvalp(e);
460 gel(z,2) = icopy(p);
461 gel(z,3) = icopy(mod);
462 gel(z,4) = icopy(u); return z;
463 }
464
465 /* Mod(x,X) + Mod(y,X) */
466 #define add_polmod_same add_polmod_scal
467 /* Mod(x,X) + Mod(y,Y) */
468 static GEN
add_polmod(GEN X,GEN Y,GEN x,GEN y)469 add_polmod(GEN X, GEN Y, GEN x, GEN y)
470 {
471 long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
472 GEN z = cgetg(3,t_POLMOD);
473 long vx = varn(X), vy = varn(Y);
474 if (vx==vy) {
475 pari_sp av;
476 gel(z,1) = srgcd(X,Y); av = avma;
477 gel(z,2) = gerepileupto(av, gmod(gadd(x, y), gel(z,1))); return z;
478 }
479 if (varncmp(vx, vy) < 0)
480 { gel(z,1) = gcopy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
481 else
482 { gel(z,1) = gcopy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
483 gel(z,2) = gadd(x, y); return z;
484 }
485 /* Mod(y, Y) + x, assuming x scalar or polynomial in same var and reduced degree */
486 static GEN
add_polmod_scal(GEN Y,GEN y,GEN x)487 add_polmod_scal(GEN Y, GEN y, GEN x)
488 {
489 GEN z = cgetg(3,t_POLMOD);
490 gel(z,1) = gcopy(Y);
491 gel(z,2) = gadd(x, y); return z;
492 }
493
494 /* check y[a..b-1] and set signe to 1 if one coeff is non-0, 0 otherwise
495 * For t_POL and t_SER */
496 static GEN
NORMALIZE_i(GEN y,long a,long b)497 NORMALIZE_i(GEN y, long a, long b)
498 {
499 long i;
500 for (i = a; i < b; i++)
501 if (!gcmp0(gel(y,i))) { setsigne(y, 1); return y; }
502 setsigne(y, 0); return y;
503 }
504 /* typ(y) == t_POL, varn(y) = vy, x "scalar" [e.g object in lower variable] */
505 static GEN
add_pol_scal(GEN y,GEN x,long vy)506 add_pol_scal(GEN y, GEN x, long vy)
507 {
508 long i, ly = lg(y);
509 GEN z;
510 if (ly <= 3) {
511 if (ly == 2) return isexactzero(x)? zeropol(vy): scalarpol(x,vy);
512 z = cgetg(3, t_POL); z[1] = y[1];
513 gel(z,2) = gadd(x,gel(y,2));
514 if (gcmp0(gel(z,2))) {
515 if (isexactzero(gel(z,2))) { avma = (pari_sp)(z+3); return zeropol(vy); }
516 setsigne(z, 0);
517 }
518 return z;
519 }
520 z = cgetg(ly,t_POL); z[1] = y[1];
521 gel(z,2) = gadd(x,gel(y,2));
522 for (i = 3; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
523 return NORMALIZE_i(z, 2, ly);
524 }
525 /* typ(y) == t_SER, varn(y) = vy, x "scalar" [e.g object in lower variable]
526 * l = valp(y) */
527 static GEN
add_ser_scal(GEN y,GEN x,long vy,long l)528 add_ser_scal(GEN y, GEN x, long vy, long l)
529 {
530 long i, j, ly;
531 pari_sp av;
532 GEN z;
533
534 if (isexactzero(x)) return gcopy(y);
535 ly = lg(y);
536 if (l < 3-ly) return gcopy(y);
537 if (l < 0)
538 {
539 z = cgetg(ly,t_SER); z[1] = y[1];
540 for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
541 gel(z,i) = gadd(x,gel(y,i)); i++;
542 for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
543 return NORMALIZE_i(z, 2, ly);
544 }
545 if (l > 0)
546 {
547 ly += l; y -= l; z = cgetg(ly,t_SER);
548 z[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
549 gel(z,2) = gcopy(x);
550 for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
551 for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
552 if (gcmp0(x)) return NORMALIZE_i(z, l+2, ly);
553 return z;
554 }
555 /* l = 0, !isexactzero(x) */
556 if (ly==2) return zeroser(vy, 0); /* 1 + O(1) --> O(1) */
557
558 av = avma; z = cgetg(ly,t_SER);
559 x = gadd(x, gel(y,2));
560 if (!isexactzero(x))
561 {
562 z[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
563 gel(z,2) = x;
564 for (i=3; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
565 if (gcmp0(x)) return NORMALIZE_i(z, 3, ly);
566 return z;
567 }
568 avma = av; /* first coeff is 0 */
569 i = 3; while (i < ly && isexactzero(gel(y,i))) i++;
570 i -= 2; ly -= i; y += i;
571 z = cgetg(ly,t_SER); z[1] = evalvalp(i) | evalvarn(vy);
572 for (j = 2; j < ly; j++) gel(z,j) = gcopy(gel(y,j));
573 return NORMALIZE_i(z, 2, ly);
574 }
575 /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
576 static GEN
add_rfrac_scal(GEN y,GEN x)577 add_rfrac_scal(GEN y, GEN x)
578 {
579 pari_sp av = avma;
580 GEN n = gadd(gmul(x, gel(y,2)), gel(y,1));
581 return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
582 }
583
584 /* x "scalar", ty != t_MAT and non-scalar */
585 static GEN
add_scal(GEN y,GEN x,long ty,long vy)586 add_scal(GEN y, GEN x, long ty, long vy)
587 {
588 long tx;
589 switch(ty)
590 {
591 case t_POL: return add_pol_scal(y, x, vy);
592 case t_SER: return add_ser_scal(y, x, vy, valp(y));
593 case t_RFRAC: return add_rfrac_scal(y, x);
594 case t_VEC: case t_COL:
595 tx = typ(x);
596 if (!is_matvec_t(tx) && isexactzero(x)) return gcopy(y);
597 break;
598 }
599 pari_err(operf,"+",x,y);
600 return NULL; /* not reached */
601 }
602
603 static GEN
addfrac(GEN x,GEN y)604 addfrac(GEN x, GEN y)
605 {
606 GEN x1 = gel(x,1), x2 = gel(x,2), z = cgetg(3,t_FRAC);
607 GEN y1 = gel(y,1), y2 = gel(y,2), p1, r, n, d, delta;
608
609 delta = gcdii(x2,y2);
610 if (is_pm1(delta))
611 { /* numerator is non-zero */
612 gel(z,1) = gerepileuptoint((pari_sp)z, addii(mulii(x1,y2), mulii(y1,x2)));
613 gel(z,2) = mulii(x2,y2); return z;
614 }
615 x2 = diviiexact(x2,delta);
616 y2 = diviiexact(y2,delta);
617 n = addii(mulii(x1,y2), mulii(y1,x2));
618 if (!signe(n)) { avma = (pari_sp)(z+3); return gen_0; }
619 d = mulii(x2, y2);
620 p1 = dvmdii(n, delta, &r);
621 if (r == gen_0)
622 {
623 if (is_pm1(d)) { avma = (pari_sp)(z+3); return icopy(p1); }
624 avma = (pari_sp)z;
625 gel(z,2) = icopy(d);
626 gel(z,1) = icopy(p1); return z;
627 }
628 p1 = gcdii(delta, r);
629 if (!is_pm1(p1))
630 {
631 delta = diviiexact(delta, p1);
632 n = diviiexact(n, p1);
633 }
634 d = mulii(d,delta); avma = (pari_sp)z;
635 gel(z,1) = icopy(n);
636 gel(z,2) = icopy(d); return z;
637 }
638
639 static GEN
add_rfrac(GEN x,GEN y)640 add_rfrac(GEN x, GEN y)
641 {
642 GEN x1 = gel(x,1), x2 = gel(x,2), z = cgetg(3,t_RFRAC);
643 GEN y1 = gel(y,1), y2 = gel(y,2), p1, r, n, d, delta;
644 pari_sp tetpil;
645
646 delta = ggcd(x2,y2);
647 if (gcmp1(delta))
648 { /* numerator is non-zero */
649 gel(z,1) = gerepileupto((pari_sp)z, gadd(gmul(x1,y2), gmul(y1,x2)));
650 gel(z,2) = gmul(x2, y2); return z;
651 }
652 x2 = gdeuc(x2,delta);
653 y2 = gdeuc(y2,delta);
654 n = gadd(gmul(x1,y2), gmul(y1,x2));
655 if (gcmp0(n)) return gerepileupto((pari_sp)(z+3), n);
656 tetpil = avma; d = gmul(x2, y2);
657 p1 = poldivrem(n, delta, &r); /* we want gcd(n,delta) */
658 if (gcmp0(r))
659 {
660 if (lg(d) == 3) /* "constant" denominator */
661 {
662 d = gel(d,2);
663 if (gcmp_1(d)) p1 = gneg(p1);
664 else if (!gcmp1(d)) p1 = gdiv(p1, d);
665 return gerepileupto((pari_sp)(z+3), p1);
666 }
667 gel(z,1) = p1; gel(z,2) = d;
668 gerepilecoeffssp((pari_sp)z,tetpil,z+1,2); return z;
669 }
670 p1 = ggcd(delta, r);
671 if (gcmp1(p1))
672 {
673 tetpil = avma;
674 gel(z,1) = gcopy(n);
675 }
676 else
677 {
678 delta = gdeuc(delta, p1);
679 tetpil = avma;
680 gel(z,1) = gdeuc(n,p1);
681 }
682 gel(z,2) = gmul(d,delta);
683 gerepilecoeffssp((pari_sp)z,tetpil,z+1,2); return z;
684 }
685
686 GEN
gadd(GEN x,GEN y)687 gadd(GEN x, GEN y)
688 {
689 long tx = typ(x), ty = typ(y), vx, vy, lx, ly, i, l;
690 pari_sp av, tetpil;
691 GEN z, p1;
692
693 if (tx == ty) switch(tx) /* shortcut to generic case */
694 {
695 case t_INT: return addii(x,y);
696 case t_REAL: return addrr(x,y);
697 case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
698 z = cgetg(3,t_INTMOD);
699 if (X==Y || equalii(X,Y))
700 return add_intmod_same(z, X, gel(x,2), gel(y,2));
701 gel(z,1) = gcdii(X,Y);
702 av = avma; p1 = addii(gel(x,2),gel(y,2));
703 gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
704 }
705 case t_FRAC: return addfrac(x,y);
706 case t_COMPLEX: z = cgetg(3,t_COMPLEX);
707 gel(z,2) = gadd(gel(x,2),gel(y,2));
708 if (isexactzero(gel(z,2)))
709 {
710 avma = (pari_sp)(z+3);
711 return gadd(gel(x,1),gel(y,1));
712 }
713 gel(z,1) = gadd(gel(x,1),gel(y,1));
714 return z;
715 case t_PADIC:
716 if (!equalii(gel(x,2),gel(y,2))) pari_err(operi,"+",x,y);
717 return addpp(x,y);
718 case t_QUAD: z = cgetg(4,t_QUAD);
719 if (!gequal(gel(x,1),gel(y,1))) pari_err(operi,"+",x,y);
720 gel(z,1) = gcopy(gel(x,1));
721 gel(z,2) = gadd(gel(x,2),gel(y,2));
722 gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
723 case t_POLMOD:
724 if (gequal(gel(x,1), gel(y,1)))
725 return add_polmod_same(gel(x,1), gel(x,2), gel(y,2));
726 return add_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
727 case t_POL:
728 vx = varn(x);
729 vy = varn(y);
730 if (vx != vy) {
731 if (varncmp(vx, vy) < 0) return add_pol_scal(x, y, vx);
732 else return add_pol_scal(y, x, vy);
733 }
734 /* same variable */
735 lx = lg(x);
736 ly = lg(y);
737 if (lx == ly) {
738 z = cgetg(lx, t_POL); z[1] = x[1];
739 for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
740 return normalizepol_i(z, lx);
741 }
742 if (ly < lx) {
743 z = cgetg(lx,t_POL); z[1] = x[1];
744 for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
745 for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
746 if (!signe(x)) z = normalizepol_i(z, lx);
747 } else {
748 z = cgetg(ly,t_POL); z[1] = y[1];
749 for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
750 for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
751 if (!signe(y)) z = normalizepol_i(z, ly);
752 }
753 return z;
754 case t_SER:
755 vx = varn(x);
756 vy = varn(y);
757 if (vx != vy) {
758 if (varncmp(vx, vy) < 0) return add_ser_scal(x, y, vx, valp(x));
759 else return add_ser_scal(y, x, vy, valp(y));
760 }
761 l = valp(y) - valp(x);
762 if (l < 0) { l = -l; swap(x,y); }
763 /* valp(x) <= valp(y) */
764 lx = lg(x);
765 ly = lg(y) + l; if (lx < ly) ly = lx;
766 if (l)
767 {
768 if (l > lx-2) return gcopy(x);
769 z = cgetg(ly,t_SER);
770 for (i=2; i<=l+1; i++) gel(z,i) = gcopy(gel(x,i));
771 for ( ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-l));
772 } else {
773 z = cgetg(ly,t_SER);
774 for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
775 }
776 z[1] = x[1]; return normalize(z);
777 case t_RFRAC:
778 vx = gvar(x);
779 vy = gvar(y);
780 if (vx != vy) {
781 if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
782 else return add_rfrac_scal(y, x);
783 }
784 return add_rfrac(x,y);
785 case t_VEC: case t_COL: case t_MAT:
786 ly = lg(y);
787 if (ly != lg(x)) pari_err(operi,"+",x,y);
788 z = cgetg(ly, ty);
789 for (i = 1; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
790 return z;
791
792 default: pari_err(operf,"+",x,y);
793 }
794 /* tx != ty */
795 if (tx > ty) { swap(x,y); lswap(tx,ty); }
796
797 if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
798 {
799 case t_INT:
800 switch(ty)
801 {
802 case t_REAL: return addir(x,y);
803 case t_INTMOD:
804 z = cgetg(3, t_INTMOD);
805 return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
806 case t_FRAC: z = cgetg(3,t_FRAC);
807 gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
808 gel(z,2) = icopy(gel(y,2)); return z;
809 case t_COMPLEX: return addRc(x, y);
810 case t_PADIC: return addQp(x,y);
811 case t_QUAD: return addRq(x, y);
812 }
813
814 case t_REAL:
815 switch(ty)
816 {
817 case t_FRAC:
818 if (!signe(y[1])) return rcopy(x);
819 if (!signe(x))
820 {
821 lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
822 return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
823 }
824 av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
825 return gerepile(av,tetpil,divri(z,gel(y,2)));
826 case t_COMPLEX: return addRc(x, y);
827 case t_QUAD: return gcmp0(y)? rcopy(x): addqf(y, x, lg(x));
828
829 default: pari_err(operf,"+",x,y);
830 }
831
832 case t_INTMOD:
833 switch(ty)
834 {
835 case t_FRAC: { GEN X = gel(x,1);
836 z = cgetg(3, t_INTMOD);
837 p1 = modii(mulii(gel(y,1), Fp_inv(gel(y,2),X)), X);
838 return add_intmod_same(z, X, p1, gel(x,2));
839 }
840 case t_COMPLEX: return addRc(x, y);
841 case t_PADIC: { GEN X = gel(x,1);
842 z = cgetg(3, t_INTMOD);
843 return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
844 }
845 case t_QUAD: return addRq(x, y);
846 }
847
848 case t_FRAC:
849 switch (ty)
850 {
851 case t_COMPLEX: return addRc(x, y);
852 case t_PADIC: return addQp(x,y);
853 case t_QUAD: return addRq(x, y);
854 }
855
856 case t_COMPLEX:
857 switch(ty)
858 {
859 case t_PADIC:
860 return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
861 case t_QUAD:
862 lx = precision(x); if (!lx) pari_err(operi,"+",x,y);
863 return gcmp0(y)? rcopy(x): addqf(y, x, lx);
864 }
865
866 case t_PADIC: /* ty == t_QUAD */
867 return (kro_quad(gel(y,1),gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
868 }
869 /* tx < ty, !is_const_t(y) */
870 if (ty == t_MAT) {
871 if (is_matvec_t(tx)) pari_err(operf,"+",x,y);
872 if (isexactzero(x)) return gcopy(y);
873 return gaddmat(x,y);
874 }
875 if (ty == t_POLMOD) /* is_const_t(tx) in this case */
876 return add_polmod_scal(gel(y,1), gel(y,2), x);
877 vy = gvar(y);
878 if (is_scalar_t(tx)) {
879 if (tx == t_POLMOD)
880 {
881 vx = varn(x[1]);
882 if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
883 else
884 if (varncmp(vx,vy) > 0) return add_scal(y, x, ty, vy);
885 return add_polmod_scal(gel(x,1), gel(x,2), y);
886 }
887 return add_scal(y, x, ty, vy);
888 }
889 /* x and y are not scalars, ty != t_MAT */
890 vx = gvar(x);
891 if (vx != vy) { /* x or y is treated as a scalar */
892 if (is_vec_t(tx) || is_vec_t(ty)) pari_err(operf,"+",x,y);
893 return (varncmp(vx, vy) < 0)? add_scal(x, y, tx, vx)
894 : add_scal(y, x, ty, vy);
895 }
896 /* vx = vy */
897 switch(tx)
898 {
899 case t_POL:
900 switch (ty)
901 {
902 case t_SER:
903 if (lg(x) == 2) return gcopy(y);
904 i = lg(y) + valp(y) - polvaluation(x, NULL);
905 if (i < 3) return gcopy(y);
906
907 p1 = greffe(x,i,0); y = gadd(p1,y);
908 free(p1); return y;
909
910 case t_RFRAC: return add_rfrac_scal(y, x);
911 }
912 break;
913
914 case t_SER:
915 if (ty == t_RFRAC)
916 {
917 GEN n, d;
918 long vn, vd;
919 n = gel(y,1); vn = gval(n, vy);
920 d = gel(y,2); vd = polvaluation(d, &d);
921
922 l = lg(x) + valp(x) - (vn - vd);
923 if (l < 3) return gcopy(x);
924
925 av = avma;
926 /* take advantage of y = t^n ! */
927 if (degpol(d)) {
928 y = gdiv(n, greffe(d,l,1));
929 } else {
930 y = gdiv(n, gel(d,2));
931 if (gvar(y) == vy) y = greffe(y,l,1); else y = scalarser(y, vy, l);
932 }
933 setvalp(y, valp(y) - vd);
934 return gerepileupto(av, gadd(y, x));
935 }
936 break;
937 }
938 pari_err(operf,"+",x,y);
939 return NULL; /* not reached */
940 }
941
942 GEN
gaddsg(long x,GEN y)943 gaddsg(long x, GEN y)
944 {
945 long ty = typ(y);
946 GEN z;
947
948 switch(ty)
949 {
950 case t_INT: return addsi(x,y);
951 case t_REAL: return addsr(x,y);
952 case t_INTMOD:
953 z = cgetg(3, t_INTMOD);
954 return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
955 case t_FRAC: z = cgetg(3,t_FRAC);
956 gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
957 gel(z,2) = icopy(gel(y,2)); return z;
958 case t_COMPLEX:
959 z = cgetg(3, t_COMPLEX);
960 gel(z,1) = gaddsg(x, gel(y,1));
961 gel(z,2) = gcopy(gel(y,2)); return z;
962
963 default: return gadd(stoi(x), y);
964 }
965 }
966
967 /********************************************************************/
968 /** **/
969 /** MULTIPLICATION **/
970 /** **/
971 /********************************************************************/
972 static GEN
mul_ser_scal(GEN y,GEN x)973 mul_ser_scal(GEN y, GEN x) {
974 long ly, i;
975 GEN z;
976 if (isexactzero(x)) { long vy = varn(y); return zeropol(vy); }
977 ly = lg(y); z = cgetg(ly,t_SER); z[1] = y[1];
978 for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
979 return normalize(z);
980 }
981 /* (n/d) * x, x "scalar" or polynomial in the same variable as d
982 * [n/d a valid RFRAC] */
983 static GEN
mul_rfrac_scal(GEN n,GEN d,GEN x)984 mul_rfrac_scal(GEN n, GEN d, GEN x)
985 {
986 pari_sp av = avma;
987 GEN z;
988 switch(typ(x))
989 {
990 case t_INTMOD: case t_POLMOD:
991 n = gmul(n, x);
992 d = gmul(d, gmodulo(gen_1, gel(x,1)));
993 return gerepileupto(av, gdiv(n,d));
994 }
995 z = gred_rfrac2_i(x, d);
996 n = simplify_i(n);
997 if (typ(z) == t_RFRAC)
998 z = gred_rfrac_simple(gmul(gel(z,1), n), gel(z,2));
999 else
1000 z = gmul(z, n);
1001 return gerepileupto(av, z);
1002 }
1003 static GEN
mul_scal(GEN y,GEN x,long ty)1004 mul_scal(GEN y, GEN x, long ty)
1005 {
1006 switch(ty)
1007 {
1008 case t_POL: return RgX_Rg_mul(y, x);
1009 case t_SER: return mul_ser_scal(y, x);
1010 case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1011 case t_QFI: case t_QFR:
1012 if (typ(x) == t_INT && gcmp1(x)) return gcopy(y); /* fall through */
1013 }
1014 pari_err(operf,"*",x,y);
1015 return NULL; /* not reached */
1016 }
1017
1018 static GEN
mul_gen_rfrac(GEN X,GEN Y)1019 mul_gen_rfrac(GEN X, GEN Y)
1020 {
1021 GEN y1 = gel(Y,1), y2 = gel(Y,2);
1022 long vx = gvar(X), vy = varn(y2);
1023 return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
1024 gred_rfrac_simple(gmul(y1,X), y2);
1025 }
1026 /* (x1/x2) * (y1/y2) */
1027 static GEN
mul_rfrac(GEN x1,GEN x2,GEN y1,GEN y2)1028 mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
1029 {
1030 GEN z, X, Y;
1031 pari_sp av = avma;
1032
1033 X = gred_rfrac2_i(x1, y2);
1034 Y = gred_rfrac2_i(y1, x2);
1035 if (typ(X) == t_RFRAC)
1036 {
1037 if (typ(Y) == t_RFRAC) {
1038 x1 = gel(X,1);
1039 x2 = gel(X,2);
1040 y1 = gel(Y,1);
1041 y2 = gel(Y,2);
1042 z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
1043 } else
1044 z = mul_gen_rfrac(Y, X);
1045 }
1046 else if (typ(Y) == t_RFRAC)
1047 z = mul_gen_rfrac(X, Y);
1048 else
1049 z = gmul(X, Y);
1050 return gerepileupto(av, z);
1051 }
1052 /* (x1/x2) /y2 */
1053 static GEN
div_rfrac_pol(GEN x1,GEN x2,GEN y2)1054 div_rfrac_pol(GEN x1, GEN x2, GEN y2)
1055 {
1056 GEN z, X;
1057 pari_sp av = avma;
1058
1059 X = gred_rfrac2_i(x1, y2);
1060 if (typ(X) == t_RFRAC)
1061 z = gred_rfrac_simple(gel(X,1), gmul(gel(X,2),x2));
1062 else
1063 z = mul_gen_rfrac(X, mkrfrac(gen_1, x2));
1064 return gerepileupto(av, z);
1065 }
1066
1067 /* Mod(y, Y) * x, assuming x scalar */
1068 static GEN
mul_polmod_scal(GEN Y,GEN y,GEN x)1069 mul_polmod_scal(GEN Y, GEN y, GEN x)
1070 {
1071 GEN z = cgetg(3,t_POLMOD);
1072 gel(z,1) = gcopy(Y);
1073 gel(z,2) = gmul(x,y); return z;
1074 }
1075 /* Mod(x,X) * Mod(y,X) */
1076 static GEN
mul_polmod_same(GEN X,GEN x,GEN y)1077 mul_polmod_same(GEN X, GEN x, GEN y)
1078 {
1079 GEN t, z = cgetg(3,t_POLMOD);
1080 pari_sp av;
1081 long v;
1082 gel(z,1) = gcopy(X); av = avma;
1083 t = gmul(x, y);
1084 /* gmod(t, gel(z,1))) optimised */
1085 if (typ(t) == t_POL && (v = varn(X)) == varn(t) && lg(t) >= lg(X))
1086 gel(z,2) = gerepileupto(av, RgX_divrem(t, X, ONLY_REM));
1087 else
1088 gel(z,2) = t;
1089 return z;
1090 }
1091 /* Mod(x,X) * Mod(y,Y) */
1092 static GEN
mul_polmod(GEN X,GEN Y,GEN x,GEN y)1093 mul_polmod(GEN X, GEN Y, GEN x, GEN y)
1094 {
1095 long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
1096 long vx = varn(X), vy = varn(Y);
1097 GEN z = cgetg(3,t_POLMOD);
1098 pari_sp av;
1099
1100 if (vx==vy) {
1101 gel(z,1) = srgcd(X,Y); av = avma;
1102 gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1))); return z;
1103 }
1104 if (varncmp(vx, vy) < 0)
1105 { gel(z,1) = gcopy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
1106 else
1107 { gel(z,1) = gcopy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
1108 gel(z,2) = gmul(x, y); return z;
1109 }
1110
1111 /* compatible t_VEC * t_COL, l = lg(x) = lg(y) */
1112 static GEN
VC_mul(GEN x,GEN y,long l)1113 VC_mul(GEN x, GEN y, long l)
1114 {
1115 pari_sp av = avma;
1116 GEN z = gen_0;
1117 long i;
1118 for (i=1; i<l; i++)
1119 {
1120 GEN c = gel(y,i);
1121 if (!isexactzeroscalar(c)) z = gadd(z, gmul(gel(x,i), c));
1122 }
1123 return gerepileupto(av,z);
1124 }
1125 /* compatible t_MAT * t_COL, l = lg(x) = lg(y), lz = l>1? lg(x[1]): 1 */
1126 static GEN
MC_mul(GEN x,GEN y,long l,long lz)1127 MC_mul(GEN x, GEN y, long l, long lz)
1128 {
1129 GEN z = cgetg(lz,t_COL);
1130 long i, j;
1131 for (i=1; i<lz; i++)
1132 {
1133 pari_sp av = avma;
1134 GEN t = gen_0;
1135 for (j=1; j<l; j++)
1136 {
1137 GEN c = gel(y,j);
1138 if (!isexactzeroscalar(c)) t = gadd(t, gmul(gcoeff(x,i,j), c));
1139 }
1140 gel(z,i) = gerepileupto(av,t);
1141 }
1142 return z;
1143 }
1144 /* x,y COMPLEX */
1145 static GEN
mulcc(GEN x,GEN y)1146 mulcc(GEN x, GEN y)
1147 {
1148 GEN xr = gel(x,1), xi = gel(x,2);
1149 GEN yr = gel(y,1), yi = gel(y,2);
1150 GEN p1, p2, p3, p4, z = cgetg(3,t_COMPLEX);
1151 pari_sp tetpil, av = avma;
1152 #if 1 /* 3M */
1153 p1 = gmul(xr,yr);
1154 p2 = gmul(xi,yi); p2 = gneg(p2);
1155 p3 = gmul(gadd(xr,xi), gadd(yr,yi));
1156 p4 = gadd(p2, gneg(p1));
1157 #else /* standard product */
1158 p1 = gmul(xr,yr);
1159 p2 = gmul(xi,yi); p2 = gneg(p2);
1160 p3 = gmul(xr,yi);
1161 p4 = gmul(xi,yr);
1162 #endif
1163 tetpil = avma;
1164 gel(z,1) = gadd(p1,p2);
1165 gel(z,2) = gadd(p3,p4);
1166 if (isexactzero(gel(z,2)))
1167 {
1168 cgiv(gel(z,2));
1169 return gerepileupto((pari_sp)(z+3), gel(z,1));
1170 }
1171 gerepilecoeffssp(av,tetpil, z+1,2); return z;
1172 }
1173 /* x,y PADIC */
1174 static GEN
mulpp(GEN x,GEN y)1175 mulpp(GEN x, GEN y) {
1176 long l = valp(x) + valp(y);
1177 pari_sp av;
1178 GEN z, t;
1179 if (!equalii(gel(x,2),gel(y,2))) pari_err(operi,"*",x,y);
1180 if (!signe(x[4])) return zeropadic(gel(x,2), l);
1181 if (!signe(y[4])) return zeropadic(gel(x,2), l);
1182
1183 t = (precp(x) > precp(y))? y: x;
1184 z = cgetp(t); setvalp(z,l); av = avma;
1185 affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
1186 avma = av; return z;
1187 }
1188 /* x,y QUAD */
1189 static GEN
mulqq(GEN x,GEN y)1190 mulqq(GEN x, GEN y) {
1191 GEN p1,p2,p3,p4, z = cgetg(4,t_QUAD);
1192 pari_sp av, tetpil;
1193 p1 = gel(x,1);
1194 if (!gequal(p1, gel(y,1))) pari_err(operi,"*",x,y);
1195
1196 gel(z,1) = gcopy(p1); av = avma;
1197 p2 = gmul(gel(x,2),gel(y,2));
1198 p3 = gmul(gel(x,3),gel(y,3));
1199 p4 = gmul(gneg_i(gel(p1,2)),p3);
1200
1201 if (gcmp0(gel(p1,3)))
1202 {
1203 tetpil = avma;
1204 gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
1205 av = avma;
1206 p2 = gmul(gel(x,2),gel(y,3));
1207 p3 = gmul(gel(x,3),gel(y,2)); tetpil = avma;
1208 gel(z,3) = gerepile(av,tetpil,gadd(p2,p3)); return z;
1209 }
1210
1211 p1 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
1212 tetpil = avma;
1213 gel(z,2) = gadd(p2,p4);
1214 gel(z,3) = gadd(p1,p3);
1215 gerepilecoeffssp(av,tetpil,z+2,2); return z;
1216 }
1217
1218 GEN
mulcxI(GEN x)1219 mulcxI(GEN x)
1220 {
1221 GEN z;
1222 switch(typ(x))
1223 {
1224 case t_INT: case t_REAL: case t_FRAC:
1225 return mkcomplex(gen_0, x);
1226 case t_COMPLEX:
1227 if (isexactzero(gel(x,1))) return gneg(gel(x,2));
1228 z = cgetg(3,t_COMPLEX);
1229 gel(z,1) = gneg(gel(x,2));
1230 z[2] = x[1]; return z;
1231 default:
1232 return gmul(gi, x);
1233 }
1234 }
1235 GEN
mulcxmI(GEN x)1236 mulcxmI(GEN x)
1237 {
1238 GEN z;
1239 switch(typ(x))
1240 {
1241 case t_INT: case t_REAL: case t_FRAC:
1242 return mkcomplex(gen_0, gneg(x));
1243 case t_COMPLEX:
1244 if (isexactzero(gel(x,1))) return gel(x,2);
1245 z = cgetg(3,t_COMPLEX);
1246 z[1] = x[2];
1247 gel(z,2) = gneg(gel(x,1)); return z;
1248 default:
1249 return gmul(mkcomplex(gen_0, gen_m1), x);
1250 }
1251 }
1252
1253 GEN
gmul(GEN x,GEN y)1254 gmul(GEN x, GEN y)
1255 {
1256 long tx, ty, lx, ly, vx, vy, i, j, l;
1257 pari_sp av, tetpil;
1258 GEN z, p1, p2;
1259
1260 if (x == y) return gsqr(x);
1261 tx = typ(x); ty = typ(y);
1262 if (tx == ty) switch(tx)
1263 {
1264 case t_INT: return mulii(x,y);
1265 case t_REAL: return mulrr(x,y);
1266 case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1267 z = cgetg(3,t_INTMOD);
1268 if (X==Y || equalii(X,Y))
1269 return mul_intmod_same(z, X, gel(x,2), gel(y,2));
1270 gel(z,1) = gcdii(X,Y); av = avma; p1 = mulii(gel(x,2),gel(y,2));
1271 gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
1272 }
1273 case t_FRAC:
1274 {
1275 GEN x1 = gel(x,1), x2 = gel(x,2);
1276 GEN y1 = gel(y,1), y2 = gel(y,2);
1277 z=cgetg(3,t_FRAC);
1278 p1 = gcdii(x1, y2);
1279 if (!is_pm1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
1280 p1 = gcdii(x2, y1);
1281 if (!is_pm1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
1282 tetpil = avma;
1283 gel(z,2) = mulii(x2,y2);
1284 gel(z,1) = mulii(x1,y1);
1285 fix_frac_if_int_GC(z,tetpil); return z;
1286 }
1287 case t_COMPLEX: return mulcc(x, y);
1288 case t_PADIC: return mulpp(x, y);
1289 case t_QUAD: return mulqq(x, y);
1290 case t_POLMOD:
1291 if (gequal(gel(x,1), gel(y,1)))
1292 return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
1293 return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
1294 case t_POL:
1295 vx = varn(x);
1296 vy = varn(y);
1297 if (vx != vy) {
1298 if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
1299 else return RgX_Rg_mul(y, x);
1300 }
1301 return RgX_mul(x, y);
1302
1303 case t_SER: {
1304 long mix, miy;
1305 vx = varn(x);
1306 vy = varn(y);
1307 if (vx != vy) {
1308 if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
1309 else return mul_ser_scal(y, x);
1310 }
1311 lx = lg(x); if (lx > lg(y)) { lx = lg(y); swap(x, y); }
1312 if (lx == 2) return zeroser(vx, valp(x)+valp(y));
1313 z = cgetg(lx,t_SER);
1314 z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
1315 if (lx > 200) /* threshold for 32bit coeffs: 400, 512 bits: 100 */
1316 {
1317 long ly;
1318 y = RgX_mul(ser2pol_i(x, lx), ser2pol_i(y, lx));
1319 ly = lg(y);
1320 if (ly >= lx) {
1321 for (i = 2; i < lx; i++) z[i] = y[i];
1322 } else {
1323 for (i = 2; i < ly; i++) z[i] = y[i];
1324 for ( ; i < lx; i++) gel(z,i) = gen_0;
1325 }
1326 z = normalize(z);
1327 return gerepilecopy((pari_sp)(z + lx), z);
1328 }
1329 x += 2; y += 2; z += 2; lx -= 3;
1330 p2 = (GEN)gpmalloc((lx+1)*sizeof(long));
1331 mix = miy = 0;
1332 for (i=0; i<=lx; i++)
1333 {
1334 p2[i] = !isexactzero(gel(y,i)); if (p2[i]) miy = i;
1335 if (!isexactzero(gel(x,i))) mix = i;
1336 p1 = gen_0; av = avma;
1337 for (j=i-mix; j<=min(i,miy); j++)
1338 if (p2[j]) p1 = gadd(p1, gmul(gel(y,j),gel(x,i-j)));
1339 gel(z,i) = gerepileupto(av,p1);
1340 }
1341 z -= 2; /* back to normalcy */
1342 free(p2); return normalize(z);
1343 }
1344 case t_QFI: return compimag(x,y);
1345 case t_QFR: return compreal(x,y);
1346 case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
1347 case t_MAT:
1348 ly = lg(y); if (ly == 1) return cgetg(1,t_MAT);
1349 lx = lg(x);
1350 if (lx != lg(y[1])) pari_err(operi,"*",x,y);
1351 z = cgetg(ly,t_MAT);
1352 l = (lx == 1)? 1: lg(x[1]);
1353 for (j=1; j<ly; j++) gel(z,j) = MC_mul(x, gel(y,j), lx, l);
1354 return z;
1355
1356 case t_VECSMALL: /* multiply as permutation. cf perm_mul */
1357 l = lg(x); z = cgetg(l, t_VECSMALL);
1358 if (l != lg(y)) break;
1359 for (i=1; i<l; i++)
1360 {
1361 long yi = y[i];
1362 if (yi < 1 || yi >= l) pari_err(operf,"*",x,y);
1363 z[i] = x[yi];
1364 }
1365 return z;
1366
1367
1368 default:
1369 pari_err(operf,"*",x,y);
1370 }
1371 /* tx != ty */
1372 if (is_const_t(ty) && is_const_t(tx)) {
1373 if (tx > ty) { swap(x,y); lswap(tx,ty); }
1374 switch(tx) {
1375 case t_INT:
1376 switch(ty)
1377 {
1378 case t_REAL: return mulir(x,y);
1379 case t_INTMOD:
1380 z = cgetg(3, t_INTMOD);
1381 return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
1382 case t_FRAC:
1383 if (!signe(x)) return gen_0;
1384 z=cgetg(3,t_FRAC);
1385 p1 = gcdii(x,gel(y,2));
1386 if (is_pm1(p1))
1387 {
1388 avma = (pari_sp)z;
1389 gel(z,2) = icopy(gel(y,2));
1390 gel(z,1) = mulii(gel(y,1), x);
1391 }
1392 else
1393 {
1394 x = diviiexact(x,p1); tetpil = avma;
1395 gel(z,2) = diviiexact(gel(y,2), p1);
1396 gel(z,1) = mulii(gel(y,1), x);
1397 fix_frac_if_int_GC(z,tetpil);
1398 }
1399 return z;
1400 case t_COMPLEX: return mulRc(x, y);
1401 case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
1402 case t_QUAD: return mulRq(x,y);
1403 }
1404
1405 case t_REAL:
1406 switch(ty)
1407 {
1408 case t_FRAC:
1409 if (!signe(y[1])) return gen_0;
1410 av = avma;
1411 return gerepileuptoleaf(av, divri(mulri(x,gel(y,1)), gel(y,2)));
1412 case t_COMPLEX: return mulRc(x, y);
1413 case t_QUAD: return mulqf(y, x, lg(x));
1414 default: pari_err(operf,"*",x,y);
1415 }
1416
1417 case t_INTMOD:
1418 switch(ty)
1419 {
1420 case t_FRAC: { GEN X = gel(x,1);
1421 z = cgetg(3, t_INTMOD); p1 = modii(mulii(gel(y,1), gel(x,2)), X);
1422 return div_intmod_same(z, X, p1, remii(gel(y,2), X));
1423 }
1424 case t_COMPLEX: return mulRc(x, y);
1425 case t_PADIC: { GEN X = gel(x,1);
1426 z = cgetg(3, t_INTMOD);
1427 return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
1428 }
1429 case t_QUAD: return mulRq(x, y);
1430 }
1431
1432 case t_FRAC:
1433 switch(ty)
1434 {
1435 case t_COMPLEX: return mulRc(x, y);
1436 case t_PADIC: return signe(x[1])? mulTp(x, y): gen_0;
1437 case t_QUAD: return mulRq(x, y);
1438 }
1439
1440 case t_COMPLEX:
1441 switch(ty)
1442 {
1443 case t_PADIC:
1444 return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
1445 case t_QUAD:
1446 lx = precision(x); if (!lx) pari_err(operi,"*",x,y);
1447 return mulqf(y, x, lx);
1448 }
1449
1450 case t_PADIC: /* ty == t_QUAD */
1451 return (kro_quad(gel(y,1),gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
1452 }
1453 }
1454
1455 if (is_matvec_t(ty))
1456 {
1457 ly = lg(y);
1458 if (!is_matvec_t(tx))
1459 {
1460 if (is_noncalc_t(tx)) pari_err(operf, "*",x,y); /* necessary if ly = 1 */
1461 z = cgetg(ly,ty);
1462 for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
1463 return z;
1464 }
1465 lx = lg(x);
1466
1467 switch(tx)
1468 {
1469 case t_VEC:
1470 switch(ty)
1471 {
1472 case t_COL:
1473 if (lx != ly) pari_err(operi,"*",x,y);
1474 return VC_mul(x, y, lx);
1475
1476 case t_MAT:
1477 if (ly == 1) return cgetg(1,t_VEC);
1478 if (lx != lg(y[1])) pari_err(operi,"*",x,y);
1479 z = cgetg(ly, t_VEC);
1480 for (i=1; i<ly; i++) gel(z,i) = VC_mul(x, gel(y,i), lx);
1481 return z;
1482 }
1483 break;
1484
1485 case t_COL:
1486 switch(ty)
1487 {
1488 case t_VEC:
1489 z = cgetg(ly,t_MAT);
1490 for (j=1; j<ly; j++)
1491 {
1492 GEN c = cgetg(lx, t_COL); gel(z,j) = c;
1493 for (i=1; i<lx; i++) gel(c,i) = gmul(gel(x,i), gel(y,j));
1494 }
1495 return z;
1496
1497 case t_MAT:
1498 if (ly != 1 && lg(y[1]) != 2) pari_err(operi,"*",x,y);
1499 z = cgetg(ly,t_MAT);
1500 for (i=1; i<ly; i++) gel(z,i) = gmul(gcoeff(y,1,i),x);
1501 return z;
1502 }
1503 break;
1504
1505 case t_MAT:
1506 switch(ty)
1507 {
1508 case t_VEC:
1509 if (lx != 2) pari_err(operi,"*",x,y);
1510 return gmul(gel(x,1), y);
1511
1512 case t_COL:
1513 if (lx != ly) pari_err(operi,"*",x,y);
1514 return MC_mul(x, y, lx, (lx == 1)? 1: lg(x[1]));
1515 }
1516 }
1517 }
1518 if (is_matvec_t(tx))
1519 {
1520 if (is_noncalc_t(ty)) pari_err(operf, "*",x,y); /* necessary if lx = 1 */
1521 lx = lg(x); z = cgetg(lx,tx);
1522 for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
1523 return z;
1524 }
1525 if (tx > ty) { swap(x,y); lswap(tx,ty); }
1526 /* tx < ty, !ismatvec(x and y) */
1527
1528 if (ty == t_POLMOD) /* is_const_t(tx) in this case */
1529 return mul_polmod_scal(gel(y,1), gel(y,2), x);
1530 if is_scalar_t(tx) {
1531 if (tx == t_POLMOD) {
1532 vx = varn(x[1]);
1533 vy = gvar(y);
1534 if (vx != vy) {
1535 if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
1536 return mul_polmod_scal(gel(x,1), gel(x,2), y);
1537 }
1538 /* error if ty == t_SER */
1539 av = avma; y = gmod(y, gel(x,1));
1540 return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
1541 }
1542 return mul_scal(y, x, ty);
1543 }
1544
1545 /* x and y are not scalars, nor matvec */
1546 vx = gvar(x);
1547 vy = gvar(y);
1548 if (vx != vy) { /* x or y is treated as a scalar */
1549 return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
1550 : mul_scal(y, x, ty);
1551 }
1552 /* vx = vy */
1553 switch(tx)
1554 {
1555 case t_POL:
1556 switch (ty)
1557 {
1558 case t_SER:
1559 {
1560 long vn;
1561 if (lg(x) == 2) return zeropol(vx);
1562 if (lg(y) == 2) return zeroser(vx, valp(y)+polvaluation(x,NULL));
1563 av = avma;
1564 vn = polvaluation(x, &x);
1565 avma = av;
1566 /* take advantage of x = t^n ! */
1567 if (degpol(x)) {
1568 p1 = greffe(x,lg(y),0);
1569 p2 = gmul(p1,y); free(p1);
1570 } else
1571 p2 = mul_ser_scal(y, gel(x,2));
1572 setvalp(p2, valp(p2) + vn);
1573 return p2;
1574 }
1575
1576 case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
1577 }
1578 break;
1579
1580 case t_SER:
1581 switch (ty)
1582 {
1583 case t_RFRAC:
1584 av = avma;
1585 return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
1586 }
1587 break;
1588 }
1589 pari_err(operf,"*",x,y);
1590 return NULL; /* not reached */
1591 }
1592
1593 int
ff_poltype(GEN * x,GEN * p,GEN * pol)1594 ff_poltype(GEN *x, GEN *p, GEN *pol)
1595 {
1596 GEN Q, P = *x, pr,p1,p2,y;
1597 long i, lx;
1598
1599 if (!signe(P)) return 0;
1600 lx = lg(P); Q = *pol;
1601 for (i=2; i<lx; i++)
1602 {
1603 p1 = gel(P,i); if (typ(p1) != t_POLMOD) {Q=NULL;break;}
1604 p2 = gel(p1,1);
1605 if (Q==NULL) { Q = p2; if (degpol(Q) <= 0) return 0; }
1606 else if (p2 != Q)
1607 {
1608 if (!gequal(p2, Q))
1609 {
1610 if (DEBUGMEM) pari_warn(warner,"different modulus in ff_poltype");
1611 return 0;
1612 }
1613 if (DEBUGMEM > 2) pari_warn(warner,"different pointers in ff_poltype");
1614 }
1615 }
1616 if (Q) {
1617 *x = P = to_Kronecker(P, Q);
1618 *pol = Q; lx = lg(P);
1619 }
1620 pr = *p; y = cgetg(lx, t_POL);
1621 for (i=lx-1; i>1; i--)
1622 {
1623 p1 = gel(P,i);
1624 switch(typ(p1))
1625 {
1626 case t_INTMOD: break;
1627 case t_INT:
1628 if (*p) p1 = modii(p1, *p);
1629 gel(y,i) = p1; continue;
1630 default:
1631 return (Q && !pr)? 1: 0;
1632 }
1633 p2 = gel(p1,1);
1634 if (pr==NULL) pr = p2;
1635 else if (p2 != pr)
1636 {
1637 if (!equalii(p2, pr))
1638 {
1639 if (DEBUGMEM) pari_warn(warner,"different modulus in ff_poltype");
1640 return 0;
1641 }
1642 if (DEBUGMEM > 2) pari_warn(warner,"different pointers in ff_poltype");
1643 }
1644 y[i] = p1[2];
1645 }
1646 y[1] = P[1];
1647 *x = y; *p = pr; return (Q || pr);
1648 }
1649
1650 GEN
gsqr(GEN x)1651 gsqr(GEN x)
1652 {
1653 long tx=typ(x), lx, i, j, l;
1654 pari_sp av, tetpil;
1655 GEN z, p1, p2, p3, p4;
1656
1657 if (is_scalar_t(tx))
1658 switch(tx)
1659 {
1660 case t_INT: return sqri(x);
1661 case t_REAL: return mulrr(x,x);
1662 case t_INTMOD: { GEN X = gel(x,1);
1663 z = cgetg(3,t_INTMOD);
1664 gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
1665 gel(z,1) = icopy(X); return z;
1666 }
1667 case t_FRAC: z=cgetg(3,t_FRAC);
1668 gel(z,1) = sqri(gel(x,1));
1669 gel(z,2) = sqri(gel(x,2)); return z;
1670
1671 case t_COMPLEX:
1672 if (isexactzero(gel(x,1))) {
1673 av = avma;
1674 return gerepileupto(av, gneg(gsqr(gel(x,2))));
1675 }
1676 z = cgetg(3,t_COMPLEX); av = avma;
1677 p1 = gadd(gel(x,1),gel(x,2));
1678 p2 = gadd(gel(x,1), gneg_i(gel(x,2)));
1679 p3 = gmul(gel(x,1),gel(x,2));
1680 tetpil = avma;
1681 gel(z,1) = gmul(p1,p2);
1682 gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
1683
1684 case t_PADIC:
1685 z = cgetg(5,t_PADIC);
1686 i = (equaliu(gel(x,2), 2) && signe(x[4]))? 1: 0;
1687 if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
1688 z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x) << 1);
1689 gel(z,2) = icopy(gel(x,2));
1690 gel(z,3) = shifti(gel(x,3), i); av = avma;
1691 gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
1692 return z;
1693
1694 case t_QUAD: z = cgetg(4,t_QUAD);
1695 p1 = gel(x,1);
1696 gel(z,1) = gcopy(p1); av = avma;
1697 p2 = gsqr(gel(x,2));
1698 p3 = gsqr(gel(x,3));
1699 p4 = gmul(gneg_i(gel(p1,2)),p3);
1700
1701 if (gcmp0(gel(p1,3)))
1702 {
1703 tetpil = avma;
1704 gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
1705 av = avma;
1706 p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
1707 gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
1708 }
1709
1710 p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
1711 tetpil = avma;
1712 gel(z,2) = gadd(p2,p4);
1713 gel(z,3) = gadd(p1,p3);
1714 gerepilecoeffssp(av,tetpil,z+2,2); return z;
1715
1716 case t_POLMOD:
1717 z=cgetg(3,t_POLMOD); gel(z,1) = gcopy(gel(x,1));
1718 av=avma; p1=gsqr(gel(x,2)); tetpil=avma;
1719 gel(z,2) = gerepile(av,tetpil, grem(p1,gel(z,1)));
1720 return z;
1721 }
1722
1723 switch(tx)
1724 {
1725 case t_POL:
1726 {
1727 GEN a = x, p = NULL, pol = NULL;
1728 long vx = varn(x);
1729 av = avma;
1730 if (ff_poltype(&x,&p,&pol))
1731 {
1732 z = FpX_sqr(x, p);
1733 if (p) z = FpX_to_mod(z,p);
1734 if (pol) z = from_Kronecker(z,pol);
1735 z = gerepileupto(av, z);
1736 }
1737 else { avma = av; z = RgX_sqr(a); }
1738 setvarn(z, vx); return z;
1739 }
1740
1741 case t_SER:
1742 {
1743 long mi;
1744 lx = lg(x);
1745 if (lx == 2) return zeroser(varn(x), 2*valp(x));
1746 z = cgetg(lx, t_SER);
1747 z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
1748 x += 2; z += 2; lx -= 3;
1749 p2 = (GEN)gpmalloc((lx+1)*sizeof(long));
1750 mi = 0;
1751 for (i=0; i<=lx; i++)
1752 {
1753 p2[i] = !isexactzero(gel(x,i)); if (p2[i]) mi = i;
1754 p1=gen_0; av=avma; l=((i+1)>>1) - 1;
1755 for (j=i-mi; j<=min(l,mi); j++)
1756 if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
1757 p1 = gshift(p1,1);
1758 if ((i&1) == 0 && p2[i>>1])
1759 p1 = gadd(p1, gsqr((GEN)x[i>>1]));
1760 gel(z,i) = gerepileupto(av,p1);
1761 }
1762 z -= 2; free(p2); return normalize(z);
1763 }
1764 case t_RFRAC: z = cgetg(3,t_RFRAC);
1765 gel(z,1) = gsqr(gel(x,1));
1766 gel(z,2) = gsqr(gel(x,2)); return z;
1767
1768 case t_MAT:
1769 lx = lg(x);
1770 if (lx!=1 && lx != lg(x[1])) pari_err(operi,"*",x,x);
1771 z = cgetg(lx, t_MAT);
1772 for (j=1; j<lx; j++) gel(z,j) = MC_mul(x, gel(x,j), lx, lx);
1773 return z;
1774
1775 case t_QFR: return sqcompreal(x);
1776 case t_QFI: return sqcompimag(x);
1777 case t_VECSMALL:
1778 l = lg(x); z = cgetg(l, t_VECSMALL);
1779 for (i=1; i<l; i++)
1780 {
1781 long xi = x[i];
1782 if (xi < 1 || xi >= l) pari_err(operf,"*",x,x);
1783 z[i] = x[xi];
1784 }
1785 return z;
1786 }
1787 pari_err(operf,"*",x,x);
1788 return NULL; /* not reached */
1789 }
1790
1791 /********************************************************************/
1792 /** **/
1793 /** DIVISION **/
1794 /** **/
1795 /********************************************************************/
1796 static GEN
div_rfrac_scal(GEN x,GEN y)1797 div_rfrac_scal(GEN x, GEN y)
1798 {
1799 pari_sp av = avma;
1800 return gerepileupto(av, gred_rfrac_simple(gel(x,1),gmul(y, gel(x,2))));
1801 }
1802 static GEN
div_scal_rfrac(GEN x,GEN y)1803 div_scal_rfrac(GEN x, GEN y)
1804 {
1805 GEN y1 = gel(y,1), y2 = gel(y,2);
1806 pari_sp av = avma;
1807 if (typ(y1) == t_POL && varn(y2) == varn(y1) && degpol(y1) > 0)
1808 return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
1809 return RgX_Rg_mul(y2, gdiv(x,y1));
1810 }
1811 static GEN
div_rfrac(GEN x,GEN y)1812 div_rfrac(GEN x, GEN y)
1813 { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
1814
1815 static GEN
div_ser_scal(GEN x,GEN y)1816 div_ser_scal(GEN x, GEN y) {
1817 long i, lx = lg(x);
1818 GEN z = cgetg_copy(lx, x); z[1] = x[1];
1819 for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1820 return normalize(z);
1821 }
1822 static GEN
div_T_scal(GEN x,GEN y,long tx)1823 div_T_scal(GEN x, GEN y, long tx) {
1824 switch(tx)
1825 {
1826 case t_POL: return RgX_Rg_div(x, y);
1827 case t_SER: return div_ser_scal(x, y);
1828 case t_RFRAC: return div_rfrac_scal(x,y);
1829 }
1830 pari_err(operf,"/",x,y);
1831 return NULL; /* not reached */
1832 }
1833
1834 static GEN
div_scal_pol(GEN x,GEN y)1835 div_scal_pol(GEN x, GEN y) {
1836 long ly = lg(y);
1837 pari_sp av;
1838 if (ly == 3) return gdiv(x,gel(y,2));
1839 if (isexactzero(x)) return zeropol(varn(y));
1840 av = avma;
1841 return gerepileupto(av, gred_rfrac_simple(x,y));
1842 }
1843 static GEN
div_scal_ser(GEN x,GEN y)1844 div_scal_ser(GEN x, GEN y) { /* TODO: improve */
1845 GEN z;
1846 long ly, i;
1847 if (gcmp0(x)) { pari_sp av=avma; return gerepileupto(av, gmul(x, ginv(y))); }
1848 ly = lg(y); z = (GEN)gpmalloc(ly*sizeof(long));
1849 z[0] = evaltyp(t_SER) | evallg(ly);
1850 z[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(y));
1851 gel(z,2) = x; for (i=3; i<ly; i++) gel(z,i) = gen_0;
1852 y = gdiv(z,y); free(z); return y;
1853 }
1854 static GEN
div_scal_T(GEN x,GEN y,long ty)1855 div_scal_T(GEN x, GEN y, long ty) {
1856 switch(ty)
1857 {
1858 case t_POL: return div_scal_pol(x, y);
1859 case t_SER: return div_scal_ser(x, y);
1860 case t_RFRAC: return div_scal_rfrac(x, y);
1861 }
1862 pari_err(operf,"/",x,y);
1863 return NULL; /* not reached */
1864 }
1865
1866 /* assume tx = ty = t_SER, same variable vx */
1867 static GEN
div_ser(GEN x,GEN y,long vx)1868 div_ser(GEN x, GEN y, long vx)
1869 {
1870 long i, j, l = valp(x) - valp(y), lx = lg(x), ly = lg(y);
1871 GEN y_lead, p1, p2, z;
1872 pari_sp av;
1873
1874 if (!signe(y)) pari_err(gdiver);
1875 if (lx == 2) return zeroser(vx, l);
1876 y_lead = gel(y,2);
1877 if (gcmp0(y_lead)) /* normalize denominator if leading term is 0 */
1878 {
1879 pari_warn(warner,"normalizing a series with 0 leading term");
1880 for (i=3,y++; i<ly; i++,y++)
1881 {
1882 y_lead = gel(y,2); ly--; l--;
1883 if (!gcmp0(y_lead)) break;
1884 }
1885 }
1886 if (ly < lx) lx = ly;
1887 p2 = (GEN)gpmalloc(lx*sizeof(long));
1888 for (i=3; i<lx; i++)
1889 {
1890 p1 = gel(y,i);
1891 if (isexactzero(p1)) p2[i] = 0;
1892 else
1893 {
1894 av = avma; gel(p2,i) = gclone(gneg_i(p1));
1895 avma = av;
1896 }
1897 }
1898 z = cgetg(lx,t_SER);
1899 z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
1900 gel(z,2) = gdiv(gel(x,2), y_lead);
1901 for (i=3; i<lx; i++)
1902 {
1903 av = avma; p1 = gel(x,i);
1904 for (j=2; j<i; j++)
1905 {
1906 l = i-j+2;
1907 if (p2[l]) p1 = gadd(p1, gmul(gel(z,j), gel(p2,l)));
1908 }
1909 gel(z,i) = gerepileupto(av, gdiv(p1, y_lead));
1910 }
1911 for (i=3; i<lx; i++)
1912 if (p2[i]) gunclone(gel(p2,i));
1913 free(p2); return normalize(z);
1914 }
1915 /* x,y compatible PADIC */
1916 static GEN
divpp(GEN x,GEN y)1917 divpp(GEN x, GEN y) {
1918 pari_sp av;
1919 long a, b;
1920 GEN z, M;
1921
1922 if (!signe(x[4])) return zeropadic(gel(x,2), valp(x)-valp(y));
1923 a = precp(x);
1924 b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
1925 z = cgetg(5, t_PADIC);
1926 z[1] = evalprecp(b) | evalvalp(valp(x) - valp(y));
1927 gel(z,2) = icopy(gel(x,2));
1928 gel(z,3) = icopy(M); av = avma;
1929 gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
1930 return z;
1931 }
1932
1933 GEN
gdiv(GEN x,GEN y)1934 gdiv(GEN x, GEN y)
1935 {
1936 long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
1937 pari_sp av, tetpil;
1938 GEN z, p1, p2;
1939
1940 if (tx == ty) switch(tx)
1941 {
1942 case t_INT:
1943 if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
1944 if (is_pm1(x)) {
1945 long s = signe(y);
1946 if (!s) pari_err(gdiver);
1947 if (signe(x) < 0) s = -s;
1948 z = cgetg(3, t_FRAC);
1949 gel(z,1) = s<0? gen_m1: gen_1;
1950 gel(z,2) = absi(y); return z;
1951 }
1952 return gred_frac2(x,y);
1953
1954 case t_REAL: return divrr(x,y);
1955 case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
1956 z = cgetg(3,t_INTMOD);
1957 if (X==Y || equalii(X,Y))
1958 return div_intmod_same(z, X, gel(x,2), gel(y,2));
1959 gel(z,1) = gcdii(X,Y); av = avma;
1960 p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
1961 gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
1962 }
1963 case t_FRAC: {
1964 GEN x1 = gel(x,1), x2 = gel(x,2);
1965 GEN y1 = gel(y,1), y2 = gel(y,2);
1966 z = cgetg(3, t_FRAC);
1967 p1 = gcdii(x1, y1);
1968 if (!is_pm1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
1969 p1 = gcdii(x2, y2);
1970 if (!is_pm1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
1971 tetpil = avma;
1972 gel(z,2) = mulii(x2,y1);
1973 gel(z,1) = mulii(x1,y2);
1974 fix_frac(z);
1975 fix_frac_if_int_GC(z,tetpil);
1976 return z;
1977 }
1978 case t_COMPLEX:
1979 av=avma; p1 = cxnorm(y); p2 = mulcc(x, gconj(y)); tetpil = avma;
1980 return gerepile(av, tetpil, gdiv(p2,p1));
1981
1982 case t_PADIC:
1983 if (!equalii(gel(x,2),gel(y,2))) pari_err(operi,"/",x,y);
1984 return divpp(x, y);
1985
1986 case t_QUAD:
1987 if (!gequal(gel(x,1),gel(y,1))) pari_err(operi,"/",x,y);
1988 av = avma; p1 = quadnorm(y); p2 = mulqq(x, gconj(y)); tetpil = avma;
1989 return gerepile(av, tetpil, gdiv(p2,p1));
1990
1991 case t_POLMOD: av = avma;
1992 if (gequal(gel(x,1), gel(y,1)))
1993 {
1994 GEN X = gel(x,1);
1995 x = gel(x,2);
1996 y = gel(y,2);
1997 if (degpol(X) == 2) { /* optimized for speed */
1998 z = mul_polmod_same(X, x, quad_polmod_conj(y, X));
1999 return gdiv(z, quad_polmod_norm(y, X));
2000 }
2001 y = ginvmod(y, X);
2002 z = mul_polmod_same(X, x, y);
2003 } else z = gmul(x, ginv(y));
2004 return gerepileupto(av, z);
2005
2006 case t_POL:
2007 vx = varn(x);
2008 vy = varn(y);
2009 if (vx != vy) {
2010 if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
2011 else return div_scal_pol(x, y);
2012 }
2013 if (!signe(y)) pari_err(gdiver);
2014 if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
2015 if (isexactzero(x)) return zeropol(vy);
2016 return gred_rfrac2(x,y);
2017
2018 case t_SER:
2019 vx = varn(x);
2020 vy = varn(y);
2021 if (vx != vy) {
2022 if (varncmp(vx, vy) < 0) return div_ser_scal(x, y);
2023 else return div_scal_ser(x, y);
2024 }
2025 return div_ser(x, y, vx);
2026 case t_RFRAC:
2027 vx = gvar(x);
2028 vy = gvar(y);
2029 if (vx != vy) {
2030 if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
2031 else return div_scal_rfrac(x, y);
2032 }
2033 return div_rfrac(x,y);
2034
2035 case t_QFI: av = avma; return gerepileupto(av, compimag(x, ginv(y)));
2036 case t_QFR: av = avma; return gerepileupto(av, compreal(x, ginv(y)));
2037
2038 case t_MAT:
2039 av = avma;
2040 return gerepileupto(av, gmul(x, invmat(y)));
2041
2042 default: pari_err(operf,"/",x,y);
2043 }
2044
2045 if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
2046 {
2047 long s = signe(x);
2048 if (!s) {
2049 if (gcmp0(y)) pari_err(gdiver);
2050 if (ty != t_INTMOD) return gen_0;
2051 z = cgetg(3,t_INTMOD);
2052 gel(z,1) = icopy(gel(y,1));
2053 gel(z,2) = gen_0; return z;
2054 }
2055 if (is_pm1(x)) {
2056 if (s > 0) return ginv(y);
2057 av = avma; return gerepileupto(av, ginv(gneg(y)));
2058 }
2059 switch(ty)
2060 {
2061 case t_REAL: return divir(x,y);
2062 case t_INTMOD:
2063 z = cgetg(3, t_INTMOD);
2064 return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
2065 case t_FRAC:
2066 z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
2067 if (is_pm1(p1))
2068 {
2069 avma = (pari_sp)z;
2070 gel(z,2) = icopy(gel(y,1));
2071 gel(z,1) = mulii(gel(y,2), x);
2072 fix_frac(z);
2073 fix_frac_if_int(z);
2074 }
2075 else
2076 {
2077 x = diviiexact(x,p1); tetpil = avma;
2078 gel(z,2) = diviiexact(gel(y,1), p1);
2079 gel(z,1) = mulii(gel(y,2), x);
2080 fix_frac(z);
2081 fix_frac_if_int_GC(z,tetpil);
2082 }
2083 return z;
2084
2085 case t_COMPLEX: return divRc(x,y);
2086 case t_PADIC: return divTp(x, y);
2087 case t_QUAD:
2088 av = avma; p1 = quadnorm(y); p2 = mulRq(x, gconj(y)); tetpil = avma;
2089 return gerepile(av, tetpil, gdiv(p2,p1));
2090 }
2091 }
2092 if (gcmp0(y) && ty != t_MAT) pari_err(gdiver);
2093
2094 if (is_const_t(tx) && is_const_t(ty)) switch(tx)
2095 {
2096 case t_REAL:
2097 switch(ty)
2098 {
2099 case t_INT: return divri(x,y);
2100 case t_FRAC:
2101 av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
2102 return gerepileuptoleaf(av, z);
2103 case t_COMPLEX: return divRc(x, y);
2104 case t_QUAD: return divfq(x, y, lg(x));
2105 default: pari_err(operf,"/",x,y);
2106 }
2107
2108 case t_INTMOD:
2109 switch(ty)
2110 {
2111 case t_INT:
2112 z = cgetg(3, t_INTMOD);
2113 return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
2114 case t_FRAC: { GEN X = gel(x,1);
2115 z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
2116 return div_intmod_same(z, X, p1, modii(gel(y,1), X));
2117 }
2118 case t_COMPLEX:
2119 av = avma; p1 = cxnorm(y); p2 = mulRc(x, gconj(y)); tetpil = avma;
2120 return gerepile(av,tetpil, gdiv(p2,p1));
2121
2122 case t_QUAD:
2123 av = avma; p1 = quadnorm(y); p2 = gmul(x,gconj(y)); tetpil = avma;
2124 return gerepile(av,tetpil, gdiv(p2,p1));
2125
2126 case t_PADIC: { GEN X = gel(x,1);
2127 z = cgetg(3, t_INTMOD);
2128 return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
2129 }
2130 case t_REAL: pari_err(operf,"/",x,y);
2131 }
2132
2133 case t_FRAC:
2134 switch(ty)
2135 {
2136 case t_INT: z = cgetg(3, t_FRAC);
2137 p1 = gcdii(y,gel(x,1));
2138 if (is_pm1(p1))
2139 {
2140 avma = (pari_sp)z; tetpil = 0;
2141 gel(z,1) = icopy(gel(x,1));
2142 }
2143 else
2144 {
2145 y = diviiexact(y,p1); tetpil = avma;
2146 gel(z,1) = diviiexact(gel(x,1), p1);
2147 }
2148 gel(z,2) = mulii(gel(x,2),y);
2149 fix_frac(z);
2150 if (tetpil) fix_frac_if_int_GC(z,tetpil);
2151 return z;
2152
2153 case t_REAL:
2154 av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
2155 return gerepile(av, tetpil, divir(gel(x,1), p1));
2156
2157 case t_INTMOD: { GEN Y = gel(y,1);
2158 z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
2159 return div_intmod_same(z, Y, gel(x,1), p1);
2160 }
2161 case t_COMPLEX: return divRc(x, y);
2162
2163 case t_PADIC:
2164 if (!signe(x[1])) return gen_0;
2165 return divTp(x, y);
2166
2167 case t_QUAD:
2168 av=avma; p1=quadnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
2169 return gerepile(av,tetpil,gdiv(p2,p1));
2170 }
2171
2172 case t_COMPLEX:
2173 switch(ty)
2174 {
2175 case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: return divcR(x,y);
2176 case t_PADIC:
2177 return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
2178 case t_QUAD:
2179 lx = precision(x); if (!lx) pari_err(operi,"/",x,y);
2180 return divfq(x, y, lx);
2181 }
2182
2183 case t_PADIC:
2184 switch(ty)
2185 {
2186 case t_INT: case t_FRAC: { GEN p = gel(x,2);
2187 return signe(x[4])? divpT(x, y)
2188 : zeropadic(p, valp(x) - ggval(y,p));
2189 }
2190 case t_INTMOD: { GEN Y = gel(y,1);
2191 z = cgetg(3, t_INTMOD);
2192 return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
2193 }
2194 case t_COMPLEX: case t_QUAD:
2195 av=avma; p1=gmul(x,gconj(y)); p2=gnorm(y); tetpil=avma;
2196 return gerepile(av,tetpil,gdiv(p1,p2));
2197
2198 case t_REAL: pari_err(operf,"/",x,y);
2199 }
2200
2201 case t_QUAD:
2202 switch (ty)
2203 {
2204 case t_INT: case t_INTMOD: case t_FRAC:
2205 z = cgetg(4,t_QUAD);
2206 gel(z,1) = gcopy(gel(x,1));
2207 gel(z,2) = gdiv(gel(x,2), y);
2208 gel(z,3) = gdiv(gel(x,3), y); return z;
2209 case t_REAL: return divqf(x, y, lg(y));
2210 case t_PADIC: return divTp(x, y);
2211 case t_COMPLEX:
2212 ly = precision(y); if (!ly) pari_err(operi,"/",x,y);
2213 return divqf(x, y, ly);
2214 }
2215 }
2216 switch(ty) {
2217 case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
2218 return gmul(x, ginv(y)); /* missing gerepile, for speed */
2219 case t_MAT:
2220 av = avma; return gerepileupto(av, gmul(x, invmat(y)));
2221 }
2222 if (is_matvec_t(tx)) {
2223 lx = lg(x); z = cgetg_copy(lx, x);
2224 for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
2225 return z;
2226 }
2227
2228 vy = gvar(y);
2229 if (tx == t_POLMOD) { GEN X = gel(x,1);
2230 vx = varn(X);
2231 if (vx != vy) {
2232 if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
2233 z = cgetg(3,t_POLMOD);
2234 gel(z,1) = gcopy(X);
2235 gel(z,2) = gdiv(gel(x,2), y); return z;
2236 }
2237 /* y is POL, SER or RFRAC */
2238 av = avma; y = ginvmod(gmod(y,X), X);
2239 return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
2240 }
2241 /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
2242 * POLMOD (done already), hence its variable is BIGINT. If the other has
2243 * variable BIGINT, then the operation is incorrect */
2244 vx = gvar(x);
2245 if (vx != vy) { /* includes cases where one is scalar */
2246 if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
2247 else return div_scal_T(x, y, ty);
2248 }
2249 switch(tx)
2250 {
2251 case t_POL:
2252 switch(ty)
2253 {
2254 case t_SER:
2255 if (lg(y) == 2)
2256 return zeroser(vx, polvaluation(x,NULL) - valp(y));
2257 p1 = greffe(x,lg(y),0);
2258 p2 = div_ser(p1, y, vx);
2259 free(p1); return p2;
2260
2261 case t_RFRAC:
2262 {
2263 GEN y1 = gel(y,1), y2 = gel(y,2);
2264 if (typ(y1) == t_POL && varn(y1) == vx)
2265 return mul_rfrac_scal(y2, y1, x);
2266 av = avma;
2267 return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
2268 }
2269 }
2270 break;
2271
2272 case t_SER:
2273 switch(ty)
2274 {
2275 case t_POL:
2276 if (lg(x) == 2)
2277 return zeroser(vx, valp(x) - polvaluation(y,NULL));
2278 p1 = greffe(y,lg(x),0);
2279 p2 = div_ser(x, p1, vx);
2280 free(p1); return p2;
2281 case t_RFRAC:
2282 av = avma;
2283 return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
2284 }
2285 break;
2286
2287 case t_RFRAC:
2288 switch(ty)
2289 {
2290 case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
2291 case t_SER:
2292 av = avma;
2293 return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2),y)));
2294 }
2295 break;
2296 }
2297 pari_err(operf,"/",x,y);
2298 return NULL; /* not reached */
2299 }
2300
2301 /********************************************************************/
2302 /** **/
2303 /** SIMPLE MULTIPLICATION **/
2304 /** **/
2305 /********************************************************************/
2306 GEN
gmulsg(long s,GEN y)2307 gmulsg(long s, GEN y)
2308 {
2309 long ty = typ(y), ly, i;
2310 pari_sp av;
2311 GEN z;
2312
2313 switch(ty)
2314 {
2315 case t_INT: return mulsi(s,y);
2316 case t_REAL: return mulsr(s,y);
2317 case t_INTMOD: { GEN p = gel(y,1);
2318 z = cgetg(3,t_INTMOD);
2319 gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
2320 gel(z,1) = icopy(p); return z;
2321 }
2322 case t_FRAC:
2323 if (!s) return gen_0;
2324 z = cgetg(3,t_FRAC);
2325 i = cgcd(s, smodis(gel(y,2), s));
2326 if (i == 1)
2327 {
2328 gel(z,2) = icopy(gel(y,2));
2329 gel(z,1) = mulis(gel(y,1), s);
2330 }
2331 else
2332 {
2333 gel(z,2) = divis(gel(y,2), i);
2334 gel(z,1) = mulis(gel(y,1), s/i);
2335 fix_frac_if_int(z);
2336 }
2337 return z;
2338
2339 case t_COMPLEX: z = cgetg(3, t_COMPLEX);
2340 gel(z,1) = gmulsg(s,gel(y,1));
2341 gel(z,2) = gmulsg(s,gel(y,2)); return z;
2342
2343 case t_PADIC:
2344 if (!s) return gen_0;
2345 av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
2346
2347 case t_QUAD: z = cgetg(4, t_QUAD);
2348 gel(z,1) = gcopy(gel(y,1));
2349 gel(z,2) = gmulsg(s,gel(y,2));
2350 gel(z,3) = gmulsg(s,gel(y,3)); return z;
2351
2352 case t_POLMOD: z = cgetg(3, t_POLMOD);
2353 gel(z,1) = gcopy(gel(y,1));
2354 gel(z,2) = gmulsg(s,gel(y,2)); return z;
2355
2356 case t_POL:
2357 if (!s || !signe(y)) return zeropol(varn(y));
2358 ly = lg(y); z = cgetg(ly,t_POL); z[1]=y[1];
2359 for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2360 return normalizepol_i(z, ly);
2361
2362 case t_SER:
2363 if (!s) return zeropol(varn(y));
2364 ly = lg(y); z = cgetg(ly,t_SER); z[1] = y[1];
2365 for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2366 return normalize(z);
2367
2368 case t_RFRAC:
2369 if (!s) return zeropol(gvar(y));
2370 z = cgetg(3, t_RFRAC);
2371 i = itos( ggcd(stoi(s),gel(y,2)) );
2372 avma = (pari_sp)z;
2373 if (i == 1)
2374 {
2375 gel(z,1) = gmulgs(gel(y,1), s);
2376 gel(z,2) = gcopy(gel(y,2));
2377 }
2378 else
2379 {
2380 gel(z,1) = gmulgs(gel(y,1), s/i);
2381 gel(z,2) = gdivgs(gel(y,2), i);
2382 }
2383 return z;
2384
2385 case t_VEC: case t_COL: case t_MAT:
2386 ly = lg(y); z = cgetg(ly,ty);
2387 for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
2388 return z;
2389 }
2390 pari_err(typeer,"gmulsg");
2391 return NULL; /* not reached */
2392 }
2393
2394 /********************************************************************/
2395 /** **/
2396 /** SIMPLE DIVISION **/
2397 /** **/
2398 /********************************************************************/
2399
2400 GEN
gdivgs(GEN x,long s)2401 gdivgs(GEN x, long s)
2402 {
2403 long tx = typ(x), lx, i;
2404 pari_sp av;
2405 GEN z, y, p1;
2406
2407 if (!s) pari_err(gdiver);
2408 switch(tx)
2409 {
2410 case t_INT:
2411 av = avma; z = divis_rem(x,s,&i);
2412 if (!i) return z;
2413
2414 i = cgcd(s, i);
2415 avma=av; z = cgetg(3,t_FRAC);
2416 if (i == 1)
2417 y = icopy(x);
2418 else
2419 {
2420 s /= i; y = diviuexact(x, i);
2421 if (signe(x) < 0) setsigne(y, -1);
2422 }
2423 gel(z,1) = y;
2424 gel(z,2) = stoi(s); fix_frac(z); return z;
2425
2426 case t_REAL:
2427 return divrs(x,s);
2428
2429 case t_INTMOD:
2430 z = cgetg(3, t_INTMOD);
2431 return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
2432
2433 case t_FRAC: z = cgetg(3, t_FRAC);
2434 i = cgcd(s, smodis(gel(x,1), s));
2435 if (i == 1)
2436 {
2437 gel(z,2) = mulsi(s, gel(x,2));
2438 gel(z,1) = icopy(gel(x,1));
2439 }
2440 else
2441 {
2442 gel(z,2) = mulsi(s/i, gel(x,2));
2443 gel(z,1) = divis(gel(x,1), i);
2444 }
2445 fix_frac(z);
2446 fix_frac_if_int(z); return z;
2447
2448 case t_COMPLEX: z = cgetg(3, t_COMPLEX);
2449 gel(z,1) = gdivgs(gel(x,1),s);
2450 gel(z,2) = gdivgs(gel(x,2),s); return z;
2451
2452 case t_PADIC: return gdiv(x, stoi(s));
2453
2454 case t_QUAD: z = cgetg(4, t_QUAD);
2455 gel(z,1) = gcopy(gel(x,1));
2456 gel(z,2) = gdivgs(gel(x,2),s);
2457 gel(z,3) = gdivgs(gel(x,3),s); return z;
2458
2459 case t_POLMOD: z = cgetg(3, t_POLMOD);
2460 gel(z,1) = gcopy(gel(x,1));
2461 gel(z,2) = gdivgs(gel(x,2),s); return z;
2462
2463 case t_RFRAC:
2464 av = avma;
2465 p1 = ggcd(stoi(s),gel(x,1));
2466 if (typ(p1) == t_INT)
2467 {
2468 avma = av;
2469 z = cgetg(3, t_RFRAC);
2470 i = p1[2];
2471 if (i == 1)
2472 {
2473 gel(z,1) = gcopy(gel(x,1));
2474 gel(z,2) = gmulsg(s,gel(x,2));
2475 }
2476 else
2477 {
2478 gel(z,1) = gdivgs(gel(x,1), i);
2479 gel(z,2) = gmulgs(gel(x,2), s/i);
2480 }
2481 }
2482 else /* t_FRAC */
2483 {
2484 z = cgetg(3, t_RFRAC);
2485 gel(z,1) = gdiv(gel(x,1), p1);
2486 gel(z,2) = gmul(gel(x,2), gdivsg(s,p1));
2487 z = gerepilecopy(av, z);
2488 }
2489 return z;
2490
2491 case t_POL: case t_SER: case t_VEC: case t_COL: case t_MAT:
2492 z = init_gen_op(x, tx, &lx, &i);
2493 for (; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
2494 return z;
2495
2496 }
2497 pari_err(operf,"/",x, stoi(s));
2498 return NULL; /* not reached */
2499 }
2500
2501 /* True shift (exact multiplication by 2^n) */
2502 GEN
gmul2n(GEN x,long n)2503 gmul2n(GEN x, long n)
2504 {
2505 long tx = typ(x), lx, i, k, l;
2506 GEN z, a, b;
2507
2508 switch(tx)
2509 {
2510 case t_INT:
2511 if (n>=0) return shifti(x,n);
2512 if (!signe(x)) return gen_0;
2513 l = vali(x); n = -n;
2514 if (n<=l) return shifti(x,-n);
2515 z = cgetg(3,t_FRAC);
2516 gel(z,1) = shifti(x,-l);
2517 gel(z,2) = int2n(n-l); return z;
2518
2519 case t_REAL:
2520 return shiftr(x,n);
2521
2522 case t_INTMOD: b = gel(x,1); a = gel(x,2);
2523 z = cgetg(3,t_INTMOD);
2524 if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
2525 gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
2526 gel(z,1) = icopy(b); return z;
2527
2528 case t_FRAC: a = gel(x,1); b = gel(x,2);
2529 l = vali(a);
2530 k = vali(b);
2531 if (n+l >= k)
2532 {
2533 if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
2534 l = n-k; k = -k;
2535 }
2536 else
2537 {
2538 k = -(l+n); l = -l;
2539 }
2540 z = cgetg(3,t_FRAC);
2541 gel(z,1) = shifti(a,l);
2542 gel(z,2) = shifti(b,k); return z;
2543
2544 case t_QUAD: z = cgetg(4,t_QUAD);
2545 gel(z,1) = gcopy(gel(x,1));
2546 gel(z,2) = gmul2n(gel(x,2),n);
2547 gel(z,3) = gmul2n(gel(x,3),n); return z;
2548
2549 case t_POLMOD: z = cgetg(3,t_POLMOD);
2550 gel(z,1) = gcopy(gel(x,1));
2551 gel(z,2) = gmul2n(gel(x,2),n); return z;
2552
2553 case t_POL:
2554 z = init_gen_op(x, tx, &lx, &i);
2555 for (; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
2556 return normalizepol_i(z, lx); /* needed if char = 2 */
2557 case t_SER:
2558 z = init_gen_op(x, tx, &lx, &i);
2559 for (; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
2560 return normalize(z); /* needed if char = 2 */
2561 case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
2562 z = init_gen_op(x, tx, &lx, &i);
2563 for (; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
2564 return z;
2565
2566 case t_RFRAC: /* int2n wrong if n < 0 */
2567 return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
2568
2569 case t_PADIC: /* int2n wrong if n < 0 */
2570 return gmul(gmul2n(gen_1,n),x);
2571 }
2572 pari_err(typeer,"gmul2n");
2573 return NULL; /* not reached */
2574 }
2575