1 /* @file poly.cc
2 *
3 * Contains the class for representing sparse multivariate
4 * polynomials with integer coefficients
5 */
6 /* #[ License : */
7 /*
8 * Copyright (C) 1984-2017 J.A.M. Vermaseren
9 * When using this file you are requested to refer to the publication
10 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
11 * This is considered a matter of courtesy as the development was paid
12 * for by FOM the Dutch physics granting agency and we would like to
13 * be able to track its scientific use to convince FOM of its value
14 * for the community.
15 *
16 * This file is part of FORM.
17 *
18 * FORM is free software: you can redistribute it and/or modify it under the
19 * terms of the GNU General Public License as published by the Free Software
20 * Foundation, either version 3 of the License, or (at your option) any later
21 * version.
22 *
23 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
24 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
25 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
26 * details.
27 *
28 * You should have received a copy of the GNU General Public License along
29 * with FORM. If not, see <http://www.gnu.org/licenses/>.
30 */
31 /* #] License : */
32 /*
33 #[ includes :
34 */
35
36 #include "poly.h"
37 #include "polygcd.h"
38
39 #include <cctype>
40 #include <cmath>
41 #include <algorithm>
42 #include <map>
43 #include <iostream>
44
45 using namespace std;
46
47 /*
48 #] includes :
49 #[ constructor (small constant polynomial) :
50 */
51
52 // constructor for a small constant polynomial
poly(PHEAD int a,WORD modp,WORD modn)53 poly::poly (PHEAD int a, WORD modp, WORD modn):
54 size_of_terms(AM.MaxTer/(LONG)sizeof(WORD)),
55 modp(0),
56 modn(1)
57 {
58 POLY_STOREIDENTITY;
59
60 terms = TermMalloc("poly::poly(int)");
61
62 if (a == 0) {
63 terms[0] = 1; // length
64 }
65 else {
66 terms[0] = 4 + AN.poly_num_vars; // length
67 terms[1] = 3 + AN.poly_num_vars; // length
68 for (int i=0; i<AN.poly_num_vars; i++) terms[2+i] = 0; // powers
69 terms[2+AN.poly_num_vars] = ABS(a); // coefficient
70 terms[3+AN.poly_num_vars] = a>0 ? 1 : -1; // length coefficient
71 }
72
73 if (modp!=-1) setmod(modp,modn);
74 }
75
76 /*
77 #] constructor (small constant polynomial) :
78 #[ constructor (large constant polynomial) :
79 */
80
81 // constructor for a large constant polynomial
poly(PHEAD const UWORD * a,WORD na,WORD modp,WORD modn)82 poly::poly (PHEAD const UWORD *a, WORD na, WORD modp, WORD modn):
83 size_of_terms(AM.MaxTer/(LONG)sizeof(WORD)),
84 modp(0),
85 modn(1)
86 {
87 POLY_STOREIDENTITY;
88
89 terms = TermMalloc("poly::poly(UWORD*,WORD)");
90
91 // remove leading zeros
92 while (*(a+ABS(na)-1)==0)
93 na -= SGN(na);
94
95 terms[0] = 3 + AN.poly_num_vars + ABS(na); // length
96 terms[1] = terms[0] - 1; // length
97 for (int i=0; i<AN.poly_num_vars; i++) terms[2+i] = 0; // powers
98 termscopy((WORD *)a, 2+AN.poly_num_vars, ABS(na)); // coefficient
99 terms[2+AN.poly_num_vars+ABS(na)] = na; // length coefficient
100
101 if (modp!=-1) setmod(modp,modn);
102 }
103
104 /*
105 #] constructor (large constant polynomial) :
106 #[ constructor (deep copy polynomial) :
107 */
108
109 // copy constructor: makes a deep copy of a polynomial
poly(const poly & a,WORD modp,WORD modn)110 poly::poly (const poly &a, WORD modp, WORD modn):
111 size_of_terms(AM.MaxTer/(LONG)sizeof(WORD))
112 {
113 POLY_GETIDENTITY(a);
114 POLY_STOREIDENTITY;
115
116 terms = TermMalloc("poly::poly(poly&)");
117
118 *this = a;
119
120 if (modp!=-1) setmod(modp,modn);
121 }
122
123 /*
124 #] constructor (deep copy polynomial) :
125 #[ destructor :
126 */
127
128 // destructor
~poly()129 poly::~poly () {
130
131 POLY_GETIDENTITY(*this);
132
133 if (size_of_terms == AM.MaxTer/(LONG)sizeof(WORD))
134 TermFree(terms, "poly::~poly");
135 else
136 M_free(terms, "poly::~poly");
137 }
138
139 /*
140 #] destructor :
141 #[ expand_memory :
142 */
143
144 // expands the memory allocated for terms to at least twice its size
expand_memory(int i)145 void poly::expand_memory (int i) {
146
147 POLY_GETIDENTITY(*this);
148
149 LONG new_size_of_terms = MaX(2*size_of_terms, i);
150 if (new_size_of_terms > MAXPOSITIVE) {
151 MLOCK(ErrorMessageLock);
152 MesPrint ((char*)"ERROR: polynomials too large (> WORDSIZE)");
153 MUNLOCK(ErrorMessageLock);
154 Terminate(1);
155 }
156
157 WORD *newterms = (WORD *)Malloc1(new_size_of_terms * sizeof(WORD), "poly::expand_memory");
158 WCOPY(newterms, terms, size_of_terms);
159
160 if (size_of_terms == AM.MaxTer/(LONG)sizeof(WORD))
161 TermFree(terms, "poly::expand_memory");
162 else
163 M_free(terms, "poly::expand_memory");
164
165 terms = newterms;
166 size_of_terms = new_size_of_terms;
167 }
168
169 /*
170 #] expand_memory :
171 #[ setmod :
172 */
173
174 // sets the coefficient space to ZZ/p^n
setmod(WORD _modp,WORD _modn)175 void poly::setmod(WORD _modp, WORD _modn) {
176
177 POLY_GETIDENTITY(*this);
178
179 if (_modp>0 && (_modp!=modp || _modn<modn)) {
180 modp = _modp;
181 modn = _modn;
182
183 WORD nmodq=0;
184 UWORD *modq=NULL;
185 bool smallq;
186
187 if (modn == 1) {
188 modq = (UWORD *)&modp;
189 nmodq = 1;
190 smallq = true;
191 }
192 else {
193 RaisPowCached(BHEAD modp,modn,&modq,&nmodq);
194 smallq = false;
195 }
196
197 coefficients_modulo(modq,nmodq,smallq);
198 }
199 else {
200 modp = _modp;
201 modn = _modn;
202 }
203 }
204
205 /*
206 #] setmod :
207 #[ coefficients_modulo :
208 */
209
210 // reduces all coefficients of the polynomial modulo a
coefficients_modulo(UWORD * a,WORD na,bool small)211 void poly::coefficients_modulo (UWORD *a, WORD na, bool small) {
212
213 POLY_GETIDENTITY(*this);
214
215 int j=1;
216 for (int i=1, di; i<terms[0]; i+=di) {
217 di = terms[i];
218
219 if (i!=j)
220 for (int k=0; k<di; k++)
221 terms[j+k] = terms[i+k];
222
223 WORD n = terms[j+terms[j]-1];
224
225 if (ABS(n)==1 && small) {
226 // small number modulo small number
227 terms[j+1+AN.poly_num_vars] = n*((UWORD)terms[j+1+AN.poly_num_vars] % a[0]);
228 if (terms[j+1+AN.poly_num_vars] == 0)
229 n=0;
230 else {
231 if (terms[j+1+AN.poly_num_vars] > +(WORD)a[0]/2) terms[j+1+AN.poly_num_vars]-=a[0];
232 if (terms[j+1+AN.poly_num_vars] < -(WORD)a[0]/2) terms[j+1+AN.poly_num_vars]+=a[0];
233 n = SGN(terms[j+1+AN.poly_num_vars]);
234 terms[j+1+AN.poly_num_vars] = ABS(terms[j+1+AN.poly_num_vars]);
235 }
236 }
237 else
238 TakeNormalModulus((UWORD *)&terms[j+1+AN.poly_num_vars], &n, a, na, NOUNPACK);
239
240 if (n!=0) {
241 terms[j] = 2+AN.poly_num_vars+ABS(n);
242 terms[j+terms[j]-1] = n;
243 j += terms[j];
244 }
245 }
246
247 terms[0] = j;
248 }
249
250 /*
251 #] coefficients_modulo :
252 #[ to_string :
253 */
254
255 // converts an integer to a string
int_to_string(WORD x)256 const string int_to_string (WORD x) {
257 char res[41];
258 snprintf (res, 41, "%i",x);
259 return res;
260 }
261
262 // converts a polynomial to a string
to_string() const263 const string poly::to_string() const {
264
265 POLY_GETIDENTITY(*this);
266
267 string res;
268
269 int printtimes;
270 UBYTE *scoeff = (UBYTE *)NumberMalloc("poly::to_string");
271
272 if (terms[0]==1)
273 // zero
274 res = "0";
275 else {
276 for (int i=1; i<terms[0]; i+=terms[i]) {
277
278 // sign
279 WORD ncoeff = terms[i+terms[i]-1];
280 if (ncoeff < 0) {
281 ncoeff*=-1;
282 res += "-";
283 }
284 else {
285 if (i>1) res += "+";
286 }
287
288 if (ncoeff==1 && terms[i+terms[i]-1-ncoeff]==1) {
289 // coeff=1, so don't print coefficient and '*'
290 printtimes = 0;
291 }
292 else {
293 // print coefficient
294 PrtLong((UWORD*)&terms[i+terms[i]-1-ncoeff], ncoeff, scoeff);
295 res += string((char *)scoeff);
296 printtimes=1;
297 }
298
299 // print variables
300 for (int j=0; j<AN.poly_num_vars; j++) {
301 if (terms[i+1+j] > 0) {
302 if (printtimes) res += "*";
303 res += string(1,'a'+j);
304 if (terms[i+1+j] > 1) res += "^" + int_to_string(terms[i+1+j]);
305 printtimes = 1;
306 }
307 }
308
309 // iff coeff=1 and all power=0, print '1' after all
310 if (!printtimes) res += "1";
311 }
312 }
313
314 // eventual modulo
315 if (modp>0) {
316 res += " (mod ";
317 res += int_to_string(modp);
318 if (modn>1) {
319 res += "^";
320 res += int_to_string(modn);
321 }
322 res += ")";
323 }
324
325 NumberFree(scoeff,"poly::to_string");
326
327 return res;
328 }
329
330 /*
331 #] to_string :
332 #[ ostream operator :
333 */
334
335 // output stream operator
operator <<(ostream & out,const poly & a)336 ostream& operator<< (ostream &out, const poly &a) {
337 return out << a.to_string();
338 }
339
340 /*
341 #] ostream operator :
342 #[ monomial_compare :
343 */
344
345 // compares two monomials (0:equal, <0:a smaller, >0:b smaller)
monomial_compare(PHEAD const WORD * a,const WORD * b)346 int poly::monomial_compare (PHEAD const WORD *a, const WORD *b) {
347 for (int i=0; i<AN.poly_num_vars; i++)
348 if (a[i+1]!=b[i+1]) return a[i+1]-b[i+1];
349 return 0;
350 }
351
352 /*
353 #] monomial_compare :
354 #[ normalize :
355 */
356
357 // brings a polynomial to normal form
normalize()358 const poly & poly::normalize() {
359
360 POLY_GETIDENTITY(*this);
361
362 WORD nmodq=0;
363 UWORD *modq=NULL;
364
365 if (modp!=0) {
366 if (modn == 1) {
367 modq = (UWORD *)&modp;
368 nmodq = 1;
369 }
370 else {
371 RaisPowCached(BHEAD modp,modn,&modq,&nmodq);
372 }
373 }
374
375 // find and sort all monomials
376 // terms[0]/num_vars+3 is an upper bound for number of terms in a
377
378 LONG poffset = AT.pWorkPointer;
379 WantAddPointers((terms[0]/(AN.poly_num_vars+3)));
380 AT.pWorkPointer += terms[0]/(AN.poly_num_vars+3);
381 #define p (&(AT.pWorkSpace[poffset]))
382
383 // WORD **p = (WORD **)Malloc1(terms[0]/(AN.poly_num_vars+3) * sizeof(WORD*), "poly::normalize");
384
385 int nterms = 0;
386 for (int i=1; i<terms[0]; i+=terms[i])
387 p[nterms++] = terms + i;
388
389 sort(p, p + nterms, monomial_larger(BHEAD0));
390
391 WORD *tmp;
392 if (size_of_terms == AM.MaxTer/(LONG)sizeof(WORD)) {
393 tmp = (WORD *)TermMalloc("poly::normalize");
394 }
395 else {
396 tmp = (WORD *)Malloc1(size_of_terms * sizeof(WORD), "poly::normalize");
397 }
398 int j=1;
399 int prevj=0;
400 tmp[0]=0;
401 tmp[1]=0;
402
403 for (int i=0; i<nterms; i++) {
404 if (i>0 && monomial_compare(BHEAD &tmp[j], p[i])==0) {
405 // duplicate term, so add coefficients
406 WORD ncoeff = tmp[j+tmp[j]-1];
407 AddLong((UWORD *)&tmp[j+1+AN.poly_num_vars], ncoeff,
408 (UWORD *)&p[i][1+AN.poly_num_vars], p[i][p[i][0]-1],
409 (UWORD *)&tmp[j+1+AN.poly_num_vars], &ncoeff);
410
411 tmp[j+1+AN.poly_num_vars+ABS(ncoeff)] = ncoeff;
412 tmp[j] = 2+AN.poly_num_vars+ABS(ncoeff);
413 }
414 else {
415 // new term
416 prevj = j;
417 j += tmp[j];
418 WCOPY(&tmp[j],p[i],p[i][0]);
419 }
420
421 if (modp!=0) {
422 // bring coefficient to normal form mod p^n
423 WORD ntmp = tmp[j+tmp[j]-1];
424 TakeNormalModulus((UWORD *)&tmp[j+1+AN.poly_num_vars], &ntmp,
425 modq,nmodq, NOUNPACK);
426 tmp[j] = 2+AN.poly_num_vars+ABS(ntmp);
427 tmp[j+tmp[j]-1] = ntmp;
428 }
429
430 // add terms to polynomial
431 if (tmp[j+tmp[j]-1]==0) {
432 tmp[j]=0;
433 j=prevj;
434 }
435 }
436
437 j+=tmp[j];
438
439 tmp[0] = j;
440 WCOPY(terms,tmp,tmp[0]);
441
442 // M_free(p, "poly::normalize");
443
444 if (size_of_terms == AM.MaxTer/(LONG)sizeof(WORD))
445 TermFree(tmp, "poly::normalize");
446 else
447 M_free(tmp, "poly::normalize");
448
449 AT.pWorkPointer = poffset;
450 #undef p
451
452 return *this;
453 }
454
455 /*
456 #] normalize :
457 #[ last_monomial_index :
458 */
459
460 // index of the last monomial, i.e., the constant term
last_monomial_index() const461 WORD poly::last_monomial_index () const {
462 POLY_GETIDENTITY(*this);
463 return terms[0] - ABS(terms[terms[0]-1]) - AN.poly_num_vars - 2;
464 }
465
466
467 // pointer to the last monomial (the constant term)
last_monomial() const468 WORD * poly::last_monomial () const {
469 return &terms[last_monomial_index()];
470 }
471
472 /*
473 #] last_monomial_index :
474 #[ compare_degree_vector :
475 */
476
compare_degree_vector(const poly & b) const477 int poly::compare_degree_vector(const poly & b) const {
478 POLY_GETIDENTITY(*this);
479
480 // special cases if one or both are 0
481 if (terms[0] == 1 && b[0] == 1) return 0;
482 if (terms[0] == 1) return -1;
483 if (b[0] == 1) return 1;
484
485 for (int i = 0; i < AN.poly_num_vars; i++) {
486 int d1 = degree(i);
487 int d2 = b.degree(i);
488 if (d1 != d2) return d1 - d2;
489 }
490
491 return 0;
492 }
493
494 /*
495 #] compare_degree_vector :
496 #[ degree_vector :
497 */
498
degree_vector() const499 std::vector<int> poly::degree_vector() const {
500 POLY_GETIDENTITY(*this);
501
502 std::vector<int> degrees(AN.poly_num_vars);
503 for (int i = 0; i < AN.poly_num_vars; i++) {
504 degrees[i] = degree(i);
505 }
506
507 return degrees;
508 }
509
510 /*
511 #] degree_vector :
512 #[ compare_degree_vector :
513 */
514
compare_degree_vector(const vector<int> & b) const515 int poly::compare_degree_vector(const vector<int> & b) const {
516 POLY_GETIDENTITY(*this);
517
518 if (terms[0] == 1) return -1;
519
520 for (int i = 0; i < AN.poly_num_vars; i++) {
521 int d1 = degree(i);
522 if (d1 != b[i]) return d1 - b[i];
523 }
524
525 return 0;
526 }
527
528 /*
529 #] compare_degree_vector :
530 #[ add :
531 */
532
533 // addition of polynomials by merging
add(const poly & a,const poly & b,poly & c)534 void poly::add (const poly &a, const poly &b, poly &c) {
535
536 POLY_GETIDENTITY(a);
537
538 c.modp = a.modp;
539 c.modn = a.modn;
540
541 WORD nmodq=0;
542 UWORD *modq=NULL;
543
544 bool both_mod_small=false;
545
546 if (c.modp!=0) {
547 if (c.modn == 1) {
548 modq = (UWORD *)&c.modp;
549 nmodq = 1;
550 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
551 both_mod_small=true;
552 }
553 else {
554 RaisPowCached(BHEAD c.modp,c.modn,&modq,&nmodq);
555 }
556 }
557
558 int ai=1,bi=1,ci=1;
559
560 while (ai<a[0] || bi<b[0]) {
561
562 c.check_memory(ci);
563
564 int cmp = ai<a[0] && bi<b[0] ? monomial_compare(BHEAD &a[ai],&b[bi]) : 0;
565
566 if (bi==b[0] || cmp>0) {
567 // insert term from a
568 c.termscopy(&a[ai],ci,a[ai]);
569 ci+=a[ai];
570 ai+=a[ai];
571 }
572 else if (ai==a[0] || cmp<0) {
573 // insert term from b
574 c.termscopy(&b[bi],ci,b[bi]);
575 ci+=b[bi];
576 bi+=b[bi];
577 }
578 else {
579 // insert term from a+b
580 c.termscopy(&a[ai],ci,MaX(a[ai],b[bi]));
581 WORD nc=0;
582 if (both_mod_small) {
583 c[ci+1+AN.poly_num_vars] = ((LONG)a[ai+1+AN.poly_num_vars]*a[ai+a[ai]-1]+
584 (LONG)b[bi+1+AN.poly_num_vars]*b[bi+b[bi]-1]) % c.modp;
585 if ((WORD)c[ci+1+AN.poly_num_vars] > +c.modp/2) c[ci+1+AN.poly_num_vars] -= c.modp;
586 if ((WORD)c[ci+1+AN.poly_num_vars] < -c.modp/2) c[ci+1+AN.poly_num_vars] += c.modp;
587 nc = (c[ci+1+AN.poly_num_vars]==0 ? 0 : SGN((WORD)c[ci+1+AN.poly_num_vars]));
588 c[ci+1+AN.poly_num_vars] = ABS((WORD)c[ci+1+AN.poly_num_vars]);
589 }
590 else {
591 AddLong((UWORD *)&a[ai+1+AN.poly_num_vars], a[ai+a[ai]-1],
592 (UWORD *)&b[bi+1+AN.poly_num_vars], b[bi+b[bi]-1],
593 (UWORD *)&c[ci+1+AN.poly_num_vars], &nc);
594 if (c.modp!=0) TakeNormalModulus((UWORD *)&c[ci+1+AN.poly_num_vars], &nc,
595 modq, nmodq, NOUNPACK);
596 }
597
598 if (nc!=0) {
599 c[ci] = 2+AN.poly_num_vars+ABS(nc);
600 c[ci+c[ci]-1] = nc;
601 ci += c[ci];
602 }
603
604 ai+=a[ai];
605 bi+=b[bi];
606 }
607 }
608
609 c[0]=ci;
610 }
611
612 /*
613 #] add :
614 #[ sub :
615 */
616
617 // subtraction of polynomials by merging
sub(const poly & a,const poly & b,poly & c)618 void poly::sub (const poly &a, const poly &b, poly &c) {
619
620 POLY_GETIDENTITY(a);
621
622 c.modp = a.modp;
623 c.modn = a.modn;
624
625 WORD nmodq=0;
626 UWORD *modq=NULL;
627
628 bool both_mod_small=false;
629
630 if (c.modp!=0) {
631 if (c.modn == 1) {
632 modq = (UWORD *)&c.modp;
633 nmodq = 1;
634 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
635 both_mod_small=true;
636 }
637 else {
638 RaisPowCached(BHEAD c.modp,c.modn,&modq,&nmodq);
639 }
640 }
641
642 int ai=1,bi=1,ci=1;
643
644 while (ai<a[0] || bi<b[0]) {
645
646 c.check_memory(ci);
647
648 int cmp = ai<a[0] && bi<b[0] ? monomial_compare(BHEAD &a[ai],&b[bi]) : 0;
649
650 if (bi==b[0] || cmp>0) {
651 // insert term from a
652 c.termscopy(&a[ai],ci,a[ai]);
653 ci+=a[ai];
654 ai+=a[ai];
655 }
656 else if (ai==a[0] || cmp<0) {
657 // insert term from b
658 c.termscopy(&b[bi],ci,b[bi]);
659 ci+=b[bi];
660 bi+=b[bi];
661 c[ci-1]*=-1;
662 }
663 else {
664 // insert term from a+b
665 c.termscopy(&a[ai],ci,MaX(a[ai],b[bi]));
666 WORD nc=0;
667 if (both_mod_small) {
668 c[ci+1+AN.poly_num_vars] = ((LONG)a[ai+1+AN.poly_num_vars]*a[ai+a[ai]-1]-
669 (LONG)b[bi+1+AN.poly_num_vars]*b[bi+b[bi]-1]) % c.modp;
670 if ((WORD)c[ci+1+AN.poly_num_vars] > +c.modp/2) c[ci+1+AN.poly_num_vars] -= c.modp;
671 if ((WORD)c[ci+1+AN.poly_num_vars] < -c.modp/2) c[ci+1+AN.poly_num_vars] += c.modp;
672 nc = (c[ci+1+AN.poly_num_vars]==0 ? 0 : SGN((WORD)c[ci+1+AN.poly_num_vars]));
673 c[ci+1+AN.poly_num_vars] = ABS((WORD)c[ci+1+AN.poly_num_vars]);
674 }
675 else {
676 AddLong((UWORD *)&a[ai+1+AN.poly_num_vars], a[ai+a[ai]-1],
677 (UWORD *)&b[bi+1+AN.poly_num_vars],-b[bi+b[bi]-1], // -b[...] causes subtraction
678 (UWORD *)&c[ci+1+AN.poly_num_vars], &nc);
679 if (c.modp!=0) TakeNormalModulus((UWORD *)&c[ci+1+AN.poly_num_vars], &nc,
680 modq, nmodq, NOUNPACK);
681 }
682
683 if (nc!=0) {
684 c[ci] = 2+AN.poly_num_vars+ABS(nc);
685 c[ci+c[ci]-1] = nc;
686 ci += c[ci];
687 }
688
689 ai+=a[ai];
690 bi+=b[bi];
691 }
692 }
693
694 c[0]=ci;
695 }
696
697 /*
698 #] sub :
699 #[ pop_heap :
700 */
701
702 // pops the largest monomial from the heap and stores it in heap[n]
pop_heap(PHEAD WORD ** heap,int n)703 void poly::pop_heap (PHEAD WORD **heap, int n) {
704
705 WORD *old = heap[0];
706
707 heap[0] = heap[--n];
708
709 int i=0;
710 while (2*i+2<n && (monomial_compare(BHEAD heap[2*i+1]+3, heap[i]+3)>0 ||
711 monomial_compare(BHEAD heap[2*i+2]+3, heap[i]+3)>0)) {
712
713 if (monomial_compare(BHEAD heap[2*i+1]+3, heap[2*i+2]+3)>0) {
714 swap(heap[i], heap[2*i+1]);
715 i=2*i+1;
716 }
717 else {
718 swap(heap[i], heap[2*i+2]);
719 i=2*i+2;
720 }
721 }
722
723 if (2*i+1<n && monomial_compare(BHEAD heap[2*i+1]+3, heap[i]+3)>0)
724 swap(heap[i], heap[2*i+1]);
725
726 heap[n] = old;
727 }
728
729 /*
730 #] pop_heap :
731 #[ push_heap :
732 */
733
734 // pushes the monomial in heap[n] onto the heap
push_heap(PHEAD WORD ** heap,int n)735 void poly::push_heap (PHEAD WORD **heap, int n) {
736
737 int i=n-1;
738
739 while (i>0 && monomial_compare(BHEAD heap[i]+3, heap[(i-1)/2]+3) > 0) {
740 swap(heap[(i-1)/2], heap[i]);
741 i=(i-1)/2;
742 }
743 }
744
745 /*
746 #] push_heap :
747 #[ mul_one_term :
748 */
749
750 // multiplies a polynomial with a monomial
mul_one_term(const poly & a,const poly & b,poly & c)751 void poly::mul_one_term (const poly &a, const poly &b, poly &c) {
752
753 POLY_GETIDENTITY(a);
754
755 int ci=1;
756
757 WORD nmodq=0;
758 UWORD *modq=NULL;
759
760 bool both_mod_small=false;
761
762 if (c.modp!=0) {
763 if (c.modn == 1) {
764 modq = (UWORD *)&c.modp;
765 nmodq = 1;
766 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
767 both_mod_small=true;
768 }
769 else {
770 RaisPowCached(BHEAD c.modp,c.modn,&modq,&nmodq);
771 }
772 }
773
774 for (int ai=1; ai<a[0]; ai+=a[ai])
775 for (int bi=1; bi<b[0]; bi+=b[bi]) {
776
777 c.check_memory(ci);
778
779 for (int i=0; i<AN.poly_num_vars; i++)
780 c[ci+1+i] = a[ai+1+i] + b[bi+1+i];
781 WORD nc;
782
783 // if both polynomials are modulo p^1, use integer calculus
784 if (both_mod_small) {
785 c[ci+1+AN.poly_num_vars] =
786 ((LONG)a[ai+1+AN.poly_num_vars] * a[ai+2+AN.poly_num_vars] *
787 b[bi+1+AN.poly_num_vars] * b[bi+2+AN.poly_num_vars]) % c.modp;
788 nc = (c[ci+1+AN.poly_num_vars]==0 ? 0 : 1);
789 }
790 else {
791 // otherwise, use form long calculus
792 MulLong((UWORD *)&a[ai+1+AN.poly_num_vars], a[ai+a[ai]-1],
793 (UWORD *)&b[bi+1+AN.poly_num_vars], b[bi+b[bi]-1],
794 (UWORD *)&c[ci+1+AN.poly_num_vars], &nc);
795 if (c.modp!=0) TakeNormalModulus((UWORD *)&c[ci+1+AN.poly_num_vars], &nc,
796 modq, nmodq, NOUNPACK);
797 }
798
799 if (nc!=0) {
800 c[ci] = 2+AN.poly_num_vars+ABS(nc);
801 ci += c[ci];
802
803 if (both_mod_small) {
804 if (c[ci-2] > +c.modp/2) c[ci-2] -= c.modp;
805 if (c[ci-2] < -c.modp/2) c[ci-2] += c.modp;
806 c[ci-1] = SGN(c[ci-2]);
807 c[ci-2] = ABS(c[ci-2]);
808 }
809 else {
810 c[ci-1] = nc;
811 }
812 }
813 }
814
815 c[0]=ci;
816 }
817
818 /*
819 #] mul_one_term :
820 #[ mul_univar :
821 */
822
823 // dense univariate multiplication, i.e., for each power find all
824 // pairs of monomials that result in that power
mul_univar(const poly & a,const poly & b,poly & c,int var)825 void poly::mul_univar (const poly &a, const poly &b, poly &c, int var) {
826
827 POLY_GETIDENTITY(a);
828
829 WORD nmodq=0;
830 UWORD *modq=NULL;
831
832 bool both_mod_small=false;
833
834 if (c.modp!=0) {
835 if (c.modn == 1) {
836 modq = (UWORD *)&c.modp;
837 nmodq = 1;
838 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
839 both_mod_small=true;
840 }
841 else {
842 RaisPowCached(BHEAD c.modp,c.modn,&modq,&nmodq);
843 }
844 }
845
846 poly t(BHEAD 0);
847 WORD nt;
848
849 int ci=1;
850
851 // bounds on the powers in a*b
852 int minpow = AN.poly_num_vars==0 ? 0 : a.last_monomial()[1+var] + b.last_monomial()[1+var];
853 int maxpow = AN.poly_num_vars==0 ? 0 : a[2+var]+b[2+var];
854 int afirst=1, blast=1;
855
856 for (int pow=maxpow; pow>=minpow; pow--) {
857
858 c.check_memory(ci);
859
860 WORD nc=0;
861
862 // adjust range in a or b
863 if (a[afirst+1+var] + b[blast+1+var] > pow) {
864 if (blast+b[blast] < b[0])
865 blast+=b[blast];
866 else
867 afirst+=a[afirst];
868 }
869
870 // find terms that result in the correct power
871 for (int ai=afirst, bi=blast; ai<a[0] && bi>=1;) {
872
873 int thispow = AN.poly_num_vars==0 ? 0 : a[ai+1+var] + b[bi+1+var];
874
875 if (thispow == pow) {
876
877 // if both polynomials are modulo p^1, use integer calculus
878 if (both_mod_small) {
879 c[ci+1+AN.poly_num_vars] =
880 ((nc==0 ? 0 : (LONG)c[ci+1+AN.poly_num_vars] * nc) +
881 (LONG)a[ai+1+AN.poly_num_vars] * a[ai+2+AN.poly_num_vars] *
882 b[bi+1+AN.poly_num_vars] * b[bi+2+AN.poly_num_vars]) % c.modp;
883 nc = (c[ci+1+AN.poly_num_vars]==0 ? 0 : 1);
884 }
885 else {
886 // otherwise, use form long calculus
887
888 MulLong((UWORD *)&a[ai+1+AN.poly_num_vars], a[ai+a[ai]-1],
889 (UWORD *)&b[bi+1+AN.poly_num_vars], b[bi+b[bi]-1],
890 (UWORD *)&t[0], &nt);
891
892 AddLong ((UWORD *)&t[0], nt,
893 (UWORD *)&c[ci+1+AN.poly_num_vars], nc,
894 (UWORD *)&c[ci+1+AN.poly_num_vars], &nc);
895
896 if (c.modp!=0) TakeNormalModulus((UWORD *)&c[ci+1+AN.poly_num_vars], &nc,
897 modq, nmodq, NOUNPACK);
898 }
899
900 ai += a[ai];
901 bi -= ABS(b[bi-1]) + 2 + AN.poly_num_vars;
902 }
903 else if (thispow > pow)
904 ai += a[ai];
905 else
906 bi -= ABS(b[bi-1]) + 2 + AN.poly_num_vars;
907 }
908
909 // add term to result
910 if (nc != 0) {
911 for (int j=0; j<AN.poly_num_vars; j++)
912 c[ci+1+j] = 0;
913 if (AN.poly_num_vars > 0)
914 c[ci+1+var] = pow;
915
916 c[ci] = 2+AN.poly_num_vars+ABS(nc);
917 ci += c[ci];
918
919 // if necessary, adjust to range [-p/2,p/2]
920 if (both_mod_small) {
921 if (c[ci-2] > +c.modp/2) c[ci-2] -= c.modp;
922 if (c[ci-2] < -c.modp/2) c[ci-2] += c.modp;
923 c[ci-1] = SGN(c[ci-2]);
924 c[ci-2] = ABS(c[ci-2]);
925 }
926 else {
927 c[ci-1] = nc;
928 }
929 }
930 }
931
932 c[0] = ci;
933 }
934
935 /*
936 #] mul_univar :
937 #[ mul_heap :
938 */
939
940 /** Multiplication of polynomials with a heap
941 *
942 * Description
943 * ===========
944 * Multiplies two multivariate polynomials. The next element of the
945 * product is efficiently determined by using a heap. If the product
946 * of the maximum power in all variables is small, a hash table is
947 * used to add equal terms for extra speed.
948 *
949 * A heap element h is formatted as follows:
950 * - h[0] = index in a
951 * - h[1] = index in b
952 * - h[2] = hash code (-1 if no hash is used)
953 * - h[3] = length of coefficient with sign
954 * - h[4...4+AN.poly_num_vars-1] = powers
955 * - h[4+AN.poly_num_vars...4+h[3]-1] = coefficient
956 */
mul_heap(const poly & a,const poly & b,poly & c)957 void poly::mul_heap (const poly &a, const poly &b, poly &c) {
958
959 POLY_GETIDENTITY(a);
960
961 WORD nmodq=0;
962 UWORD *modq=NULL;
963
964 bool both_mod_small=false;
965
966 if (c.modp!=0) {
967 if (c.modn == 1) {
968 modq = (UWORD *)&c.modp;
969 nmodq = 1;
970 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
971 both_mod_small=true;
972 }
973 else {
974 RaisPowCached(BHEAD c.modp,c.modn,&modq,&nmodq);
975 }
976 }
977
978 // find maximum powers in different variables
979 WORD *maxpower = AT.WorkPointer;
980 AT.WorkPointer += AN.poly_num_vars;
981 WORD *maxpowera = AT.WorkPointer;
982 AT.WorkPointer += AN.poly_num_vars;
983 WORD *maxpowerb = AT.WorkPointer;
984 AT.WorkPointer += AN.poly_num_vars;
985
986 for (int i=0; i<AN.poly_num_vars; i++)
987 maxpowera[i] = maxpowerb[i] = 0;
988
989 for (int ai=1; ai<a[0]; ai+=a[ai])
990 for (int j=0; j<AN.poly_num_vars; j++)
991 maxpowera[j] = MaX(maxpowera[j], a[ai+1+j]);
992
993 for (int bi=1; bi<b[0]; bi+=b[bi])
994 for (int j=0; j<AN.poly_num_vars; j++)
995 maxpowerb[j] = MaX(maxpowerb[j], b[bi+1+j]);
996
997 for (int i=0; i<AN.poly_num_vars; i++)
998 maxpower[i] = maxpowera[i] + maxpowerb[i];
999
1000 // if PROD(maxpower) small, use hash array
1001 bool use_hash = true;
1002 int nhash = 1;
1003
1004 for (int i=0; i<AN.poly_num_vars; i++) {
1005 if (nhash > POLY_MAX_HASH_SIZE / (maxpower[i]+1)) {
1006 nhash = 1;
1007 use_hash = false;
1008 break;
1009 }
1010 nhash *= maxpower[i]+1;
1011 }
1012
1013 // allocate heap and hash
1014 int nheap=a.number_of_terms();
1015
1016 WantAddPointers(nheap+nhash);
1017 WORD **heap = AT.pWorkSpace + AT.pWorkPointer;
1018
1019 for (int ai=1, i=0; ai<a[0]; ai+=a[ai], i++) {
1020 heap[i] = (WORD *) NumberMalloc("poly::mul_heap");
1021 heap[i][0] = ai;
1022 heap[i][1] = 1;
1023 heap[i][2] = -1;
1024 heap[i][3] = 0;
1025 for (int j=0; j<AN.poly_num_vars; j++)
1026 heap[i][4+j] = MAXPOSITIVE;
1027 }
1028
1029 WORD **hash = AT.pWorkSpace + AT.pWorkPointer + nheap;
1030 for (int i=0; i<nhash; i++)
1031 hash[i] = NULL;
1032
1033 int ci = 1;
1034
1035 // multiply
1036 while (nheap > 0) {
1037
1038 c.check_memory(ci);
1039
1040 pop_heap(BHEAD heap, nheap--);
1041 WORD *p = heap[nheap];
1042
1043 // if non-zero
1044 if (p[3] != 0) {
1045 if (use_hash) hash[p[2]] = NULL;
1046
1047 c[0] = ci;
1048
1049 // append this term to the result
1050 if (use_hash || ci==1 || monomial_compare(BHEAD p+3, c.last_monomial())!=0) {
1051 p[4 + AN.poly_num_vars + ABS(p[3])] = p[3];
1052 p[3] = 2 + AN.poly_num_vars + ABS(p[3]);
1053 c.termscopy(&p[3],ci,p[3]);
1054 ci += c[ci];
1055 }
1056 else {
1057 // add this term to the last term of the result
1058 ci = c.last_monomial_index();
1059 WORD nc = c[ci+c[ci]-1];
1060
1061 // if both polynomials are modulo p^1, use integer calculus
1062 if (both_mod_small) {
1063 c[ci+AN.poly_num_vars+1] = ((LONG)c[ci+AN.poly_num_vars+1]*nc + p[4+AN.poly_num_vars]*p[3]) % c.modp;
1064 if (c[ci+1+AN.poly_num_vars]==0)
1065 nc = 0;
1066 else {
1067 if (c[ci+1+AN.poly_num_vars] > +c.modp/2) c[ci+1+AN.poly_num_vars] -= c.modp;
1068 if (c[ci+1+AN.poly_num_vars] < -c.modp/2) c[ci+1+AN.poly_num_vars] += c.modp;
1069 nc = SGN(c[ci+1+AN.poly_num_vars]);
1070 c[ci+1+AN.poly_num_vars] = ABS(c[ci+1+AN.poly_num_vars]);
1071 }
1072 }
1073 else {
1074 // otherwise, use form long calculus
1075 AddLong ((UWORD *)&p[4+AN.poly_num_vars], p[3],
1076 (UWORD *)&c[ci+AN.poly_num_vars+1], nc,
1077 (UWORD *)&c[ci+AN.poly_num_vars+1],&nc);
1078
1079 if (c.modp!=0) TakeNormalModulus((UWORD *)&c[ci+1+AN.poly_num_vars], &nc,
1080 modq, nmodq, NOUNPACK);
1081 }
1082
1083 if (nc!=0) {
1084 c[ci] = 2 + AN.poly_num_vars + ABS(nc);
1085 ci += c[ci];
1086 c[ci-1] = nc;
1087 }
1088 }
1089 }
1090
1091 // add new term to the heap (ai, bi+1)
1092 while (p[1] < b[0]) {
1093
1094 for (int j=0; j<AN.poly_num_vars; j++)
1095 p[4+j] = a[p[0]+1+j] + b[p[1]+1+j];
1096
1097 // if both polynomials are modulo p^1, use integer calculus
1098 if (both_mod_small) {
1099 p[4+AN.poly_num_vars] = ((LONG)a[p[0]+1+AN.poly_num_vars]*a[p[0]+a[p[0]]-1]*
1100 b[p[1]+1+AN.poly_num_vars]*b[p[1]+b[p[1]]-1]) % c.modp;
1101 if (p[4+AN.poly_num_vars]==0)
1102 p[3]=0;
1103 else {
1104 if (p[4+AN.poly_num_vars] > +c.modp/2) p[4+AN.poly_num_vars] -= c.modp;
1105 if (p[4+AN.poly_num_vars] < -c.modp/2) p[4+AN.poly_num_vars] += c.modp;
1106 p[3] = SGN(p[4+AN.poly_num_vars]);
1107 p[4+AN.poly_num_vars] = ABS(p[4+AN.poly_num_vars]);
1108 }
1109 }
1110 else {
1111 // otherwise, use form long calculus
1112 MulLong((UWORD *)&a[p[0]+1+AN.poly_num_vars], a[p[0]+a[p[0]]-1],
1113 (UWORD *)&b[p[1]+1+AN.poly_num_vars], b[p[1]+b[p[1]]-1],
1114 (UWORD *)&p[4+AN.poly_num_vars], &p[3]);
1115
1116 if (c.modp!=0) TakeNormalModulus((UWORD *)&p[4+AN.poly_num_vars], &p[3],
1117 modq, nmodq, NOUNPACK);
1118 }
1119
1120 p[1] += b[p[1]];
1121
1122 if (use_hash) {
1123 int ID = 0;
1124 for (int i=0; i<AN.poly_num_vars; i++)
1125 ID = (maxpower[i]+1)*ID + p[4+i];
1126
1127 // if hash and unused, push onto heap
1128 if (hash[ID] == NULL) {
1129 p[2] = ID;
1130 hash[ID] = p;
1131 push_heap(BHEAD heap, ++nheap);
1132 break;
1133 }
1134 else {
1135 // if hash and used, add to heap element
1136 WORD *h = hash[ID];
1137 // if both polynomials are modulo p^1, use integer calculus
1138 if (both_mod_small) {
1139 h[4+AN.poly_num_vars] = ((LONG)p[4+AN.poly_num_vars]*p[3] + h[4+AN.poly_num_vars]*h[3]) % c.modp;
1140 if (h[4+AN.poly_num_vars]==0)
1141 h[3]=0;
1142 else {
1143 if (h[4+AN.poly_num_vars] > +c.modp/2) h[4+AN.poly_num_vars] -= c.modp;
1144 if (h[4+AN.poly_num_vars] < -c.modp/2) h[4+AN.poly_num_vars] += c.modp;
1145 h[3] = SGN(h[4+AN.poly_num_vars]);
1146 h[4+AN.poly_num_vars] = ABS(h[4+AN.poly_num_vars]);
1147 }
1148 }
1149 else {
1150 // otherwise, use form long calculus
1151 AddLong ((UWORD *)&p[4+AN.poly_num_vars], p[3],
1152 (UWORD *)&h[4+AN.poly_num_vars], h[3],
1153 (UWORD *)&h[4+AN.poly_num_vars], &h[3]);
1154
1155 if (c.modp!=0) TakeNormalModulus((UWORD *)&h[4+AN.poly_num_vars], &h[3],
1156 modq, nmodq, NOUNPACK);
1157 }
1158 }
1159 }
1160 else {
1161 // if no hash, push onto heap
1162 p[2] = -1;
1163 push_heap(BHEAD heap, ++nheap);
1164 break;
1165 }
1166 }
1167 }
1168
1169 c[0] = ci;
1170
1171 for (int ai=1, i=0; ai<a[0]; ai+=a[ai], i++)
1172 NumberFree(heap[i],"poly::mul_heap");
1173 AT.WorkPointer -= 3 * AN.poly_num_vars;
1174 }
1175
1176 /*
1177 #] mul_heap :
1178 #[ mul :
1179 */
1180
1181 /** Polynomial multiplication
1182 *
1183 * Description
1184 * ===========
1185 * This routine determines which multiplication routine to use for
1186 * multiplying two polynomials. The logic is as follows:
1187 * - If a or b consists of only one term, call mul_one_term;
1188 * - Otherwise, if both are univariate and dense, call mul_univar;
1189 * - Otherwise, call mul_heap.
1190 */
mul(const poly & a,const poly & b,poly & c)1191 void poly::mul (const poly &a, const poly &b, poly &c) {
1192
1193 c.modp = a.modp;
1194 c.modn = a.modn;
1195
1196 if (a.is_zero() || b.is_zero()) { c[0]=1; return; }
1197 if (a.is_one()) {
1198 c.check_memory(b[0]);
1199 c.termscopy(b.terms,0,b[0]);
1200 // normalize if necessary
1201 if (a.modp!=b.modp || a.modn!=b.modn) {
1202 c.modp=0;
1203 c.setmod(a.modp,a.modn);
1204 }
1205 return;
1206 }
1207 if (b.is_one()) {
1208 c.check_memory(a[0]);
1209 c.termscopy(a.terms,0,a[0]);
1210 return;
1211 }
1212
1213 int na=a.number_of_terms();
1214 int nb=b.number_of_terms();
1215
1216 if (na==1 || nb==1) {
1217 mul_one_term(a,b,c);
1218 return;
1219 }
1220
1221 int vara = a.is_dense_univariate();
1222 int varb = b.is_dense_univariate();
1223
1224 if (vara!=-2 && varb!=-2 && vara==varb) {
1225 mul_univar(a,b,c,vara);
1226 return;
1227 }
1228
1229 if (na < nb)
1230 mul_heap(a,b,c);
1231 else
1232 mul_heap(b,a,c);
1233 }
1234
1235 /*
1236 #] mul :
1237 #[ divmod_one_term : :
1238 */
1239
1240 // divides a polynomial by a monomial
divmod_one_term(const poly & a,const poly & b,poly & q,poly & r,bool only_divides)1241 void poly::divmod_one_term (const poly &a, const poly &b, poly &q, poly &r, bool only_divides) {
1242
1243 POLY_GETIDENTITY(a);
1244
1245 int qi=1, ri=1;
1246
1247 WORD nmodq=0;
1248 UWORD *modq=NULL;
1249
1250 WORD nltbinv=0;
1251 UWORD *ltbinv=NULL;
1252
1253 bool both_mod_small=false;
1254
1255 if (q.modp!=0) {
1256 if (q.modn == 1) {
1257 modq = (UWORD *)&q.modp;
1258 nmodq = 1;
1259 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
1260 both_mod_small=true;
1261 }
1262 else {
1263 RaisPowCached(BHEAD q.modp,q.modn,&modq,&nmodq);
1264 }
1265
1266 ltbinv = NumberMalloc("poly::div_one_term");
1267
1268 if (both_mod_small) {
1269 WORD ltb = b[b[1]]*b[2+AN.poly_num_vars];
1270 GetModInverses(ltb + (ltb<0?q.modp:0), q.modp, (WORD*)ltbinv, NULL);
1271 nltbinv = 1;
1272 }
1273 else
1274 GetLongModInverses(BHEAD (UWORD *)&b[2+AN.poly_num_vars], b[b[1]], modq, nmodq, ltbinv, &nltbinv, NULL, NULL);
1275 }
1276
1277 for (int ai=1; ai<a[0]; ai+=a[ai]) {
1278
1279 q.check_memory(qi);
1280 r.check_memory(ri);
1281
1282 // check divisibility of powers
1283 bool div=true;
1284 for (int j=0; j<AN.poly_num_vars; j++) {
1285 q[qi+1+j] = a[ai+1+j]-b[2+j];
1286 r[ri+1+j] = a[ai+1+j];
1287 if (q[qi+1+j] < 0) div=false;
1288 }
1289
1290 WORD nq,nr;
1291
1292 if (div) {
1293 // if variables are divisable, divide coefficient
1294 if (q.modp==0) {
1295 DivLong((UWORD *)&a[ai+1+AN.poly_num_vars], a[ai+a[ai]-1],
1296 (UWORD *)&b[2+AN.poly_num_vars], b[b[1]],
1297 (UWORD *)&q[qi+1+AN.poly_num_vars], &nq,
1298 (UWORD *)&r[ri+1+AN.poly_num_vars], &nr);
1299 }
1300 else {
1301 // if both polynomials are modulo p^1, use integer calculus
1302 if (both_mod_small) {
1303 q[qi+1+AN.poly_num_vars] =
1304 ((LONG)a[ai+1+AN.poly_num_vars] * a[ai+a[ai]-1] * ltbinv[0] * nltbinv) % q.modp;
1305 nq = (q[qi+1+AN.poly_num_vars]==0 ? 0 : 1);
1306 }
1307 else {
1308 // otherwise, use form long calculus
1309 MulLong((UWORD *)&a[ai+1+AN.poly_num_vars], a[ai+a[ai]-1],
1310 ltbinv, nltbinv,
1311 (UWORD *)&q[qi+1+AN.poly_num_vars], &nq);
1312 TakeNormalModulus((UWORD *)&q[qi+1+AN.poly_num_vars], &nq,
1313 modq,nmodq, NOUNPACK);
1314 }
1315
1316 nr=0;
1317 }
1318 }
1319 else {
1320 // if not, term becomes part of the remainder
1321 nq=0;
1322 nr=a[ai+a[ai]-1];
1323 r.termscopy(&a[ai+1+AN.poly_num_vars], ri+1+AN.poly_num_vars, ABS(nr));
1324 }
1325
1326 // add terms to quotient/remainder
1327 if (nq!=0) {
1328 q[qi] = 2+AN.poly_num_vars+ABS(nq);
1329 qi += q[qi];
1330
1331 // if necessary, adjust to range [-p/2,p/2]
1332 if (both_mod_small) {
1333 if (q[qi-2] > +q.modp/2) q[qi-2] -= q.modp;
1334 if (q[qi-2] < -q.modp/2) q[qi-2] += q.modp;
1335 q[qi-1] = SGN(q[qi-2]);
1336 q[qi-2] = ABS(q[qi-2]);
1337 }
1338 else {
1339 q[qi-1] = nq;
1340 }
1341 }
1342
1343 if (nr != 0) {
1344 if (only_divides) { r = poly(BHEAD 1); ri=r[0]; break; }
1345 r[ri] = 2+AN.poly_num_vars+ABS(nr);
1346 ri += r[ri];
1347 r[ri-1] = nr;
1348 }
1349 }
1350
1351 q[0]=qi;
1352 r[0]=ri;
1353
1354 if (q.modp!=0 || ltbinv != NULL) NumberFree(ltbinv,"poly::div_one_term");
1355 }
1356
1357 /*
1358 #] divmod_one_term :
1359 #[ divmod_univar : :
1360 */
1361
1362 /** Division of dense univariate polynomials.
1363 *
1364 * Description
1365 * ===========
1366 * Divides two dense univariate polynomials. For each power, the
1367 * method collects all terms that result in that power.
1368 *
1369 * Relevant formula [Q=A/B, P=SUM(p_i*x^i), n=deg(A), m=deg(B)]:
1370 * q_k = [ a_{m+k} - SUM(i=k+1...n-m, b_{m+k-i}*q_i) ] / b_m
1371 */
divmod_univar(const poly & a,const poly & b,poly & q,poly & r,int var,bool only_divides)1372 void poly::divmod_univar (const poly &a, const poly &b, poly &q, poly &r, int var, bool only_divides) {
1373
1374 POLY_GETIDENTITY(a);
1375
1376 WORD nmodq=0;
1377 UWORD *modq=NULL;
1378
1379 WORD nltbinv=0;
1380 UWORD *ltbinv=NULL;
1381
1382 bool both_mod_small=false;
1383
1384 if (q.modp!=0) {
1385 if (q.modn == 1) {
1386 modq = (UWORD *)&q.modp;
1387 nmodq = 1;
1388 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
1389 both_mod_small=true;
1390 }
1391 else {
1392 RaisPowCached(BHEAD q.modp,q.modn,&modq,&nmodq);
1393 }
1394 ltbinv = NumberMalloc("poly::div_univar");
1395
1396 if (both_mod_small) {
1397 WORD ltb = b[b[1]]*b[2+AN.poly_num_vars];
1398 GetModInverses(ltb + (ltb<0?q.modp:0), q.modp, (WORD*)ltbinv, NULL);
1399 nltbinv = 1;
1400 }
1401 else
1402 GetLongModInverses(BHEAD (UWORD *)&b[2+AN.poly_num_vars], b[b[1]], modq, nmodq, ltbinv, &nltbinv, NULL, NULL);
1403 }
1404
1405 WORD ns=0;
1406 WORD nt;
1407 UWORD *s = NumberMalloc("poly::div_univar");
1408 UWORD *t = NumberMalloc("poly::div_univar");
1409 s[0]=0;
1410
1411 int bpow = b[2+var];
1412
1413 int ai=1, qi=1, ri=1;
1414
1415 for (int pow=a[2+var]; pow>=0; pow--) {
1416
1417 q.check_memory(qi);
1418 r.check_memory(ri);
1419
1420 // look for the correct power in a
1421 while (ai<a[0] && a[ai+1+var] > pow)
1422 ai+=a[ai];
1423
1424 // first term of the r.h.s. of the above equation
1425 if (ai<a[0] && a[ai+1+var] == pow) {
1426 ns = a[ai+a[ai]-1];
1427 WCOPY(s, &a[ai+1+AN.poly_num_vars], ABS(ns));
1428 }
1429 else {
1430 ns = 0;
1431 }
1432
1433 int bi=1, qj=qi;
1434
1435 // second term(s) of the r.h.s. of the above equation
1436 while (qj>1 && bi<b[0]) {
1437
1438 qj -= 2 + AN.poly_num_vars + ABS(q[qj-1]);
1439
1440 while (bi<b[0] && b[bi+1+var]+q[qj+1+var] > pow)
1441 bi += b[bi];
1442
1443 if (bi<b[0] && b[bi+1+var]+q[qj+1+var] == pow) {
1444 // if both polynomials are modulo p^1, use integer calculus
1445
1446 if (both_mod_small) {
1447 s[0] = ((WORD)s[0]*ns - (LONG)b[bi+1+AN.poly_num_vars] * b[bi+b[bi]-1] *
1448 q[qj+1+AN.poly_num_vars] * q[qj+q[qj]-1]) % q.modp;
1449 ns = (s[0]==0 ? 0 : 1);
1450 }
1451 else {
1452 // otherwise, use form long calculus
1453 MulLong((UWORD *)&b[bi+1+AN.poly_num_vars], b[bi+b[bi]-1],
1454 (UWORD *)&q[qj+1+AN.poly_num_vars], q[qj+q[qj]-1],
1455 t, &nt);
1456 nt *= -1;
1457 AddLong(t,nt,s,ns,s,&ns);
1458 if (q.modp!=0) TakeNormalModulus((UWORD *)s,&ns,modq, nmodq, NOUNPACK);
1459 }
1460 }
1461 }
1462
1463 if (ns != 0) {
1464 // if necessary, adjust to range [-p/2,p/2]
1465 if (both_mod_small) {
1466 s[0]*=ns;
1467 if ((WORD)s[0] > +q.modp/2) s[0] -= q.modp;
1468 if ((WORD)s[0] < -q.modp/2) s[0] += q.modp;
1469 ns = SGN((WORD)s[0]);
1470 s[0] = ABS((WORD)s[0]);
1471 }
1472
1473 if (pow >= bpow) {
1474 // large power, so divide by b
1475 if (q.modp == 0) {
1476 DivLong(s, ns,
1477 (UWORD *)&b[2+AN.poly_num_vars], b[b[1]],
1478 (UWORD *)&q[qi+1+AN.poly_num_vars], &ns, t, &nt);
1479 }
1480 else {
1481 if (both_mod_small) {
1482 q[qi+1+AN.poly_num_vars] = ((LONG)s[0]*ns*ltbinv[0]*nltbinv) % q.modp;
1483 if ((WORD)q[qi+1+AN.poly_num_vars] > +q.modp/2) q[qi+1+AN.poly_num_vars] -= q.modp;
1484 if ((WORD)q[qi+1+AN.poly_num_vars] < -q.modp/2) q[qi+1+AN.poly_num_vars] += q.modp;
1485 ns = (q[qi+1+AN.poly_num_vars]==0 ? 0 : SGN((WORD)q[qi+1+AN.poly_num_vars]));
1486 q[qi+1+AN.poly_num_vars] = ABS((WORD)q[qi+1+AN.poly_num_vars]);
1487 }
1488 else {
1489 MulLong(s, ns, ltbinv, nltbinv, (UWORD *)&q[qi+1+AN.poly_num_vars], &ns);
1490 TakeNormalModulus((UWORD *)&q[qi+1+AN.poly_num_vars], &ns,
1491 modq,nmodq, NOUNPACK);
1492 }
1493 nt=0;
1494 }
1495 }
1496 else {
1497 // small power, so remainder
1498 WCOPY(t,s,ABS(ns));
1499 nt = ns;
1500 ns = 0;
1501 }
1502
1503 // add terms to quotient/remainder
1504 if (ns!=0) {
1505 for (int i=0; i<AN.poly_num_vars; i++)
1506 q[qi+1+i] = 0;
1507 q[qi+1+var] = pow-bpow;
1508
1509 q[qi] = 2+AN.poly_num_vars+ABS(ns);
1510 qi += q[qi];
1511 q[qi-1] = ns;
1512 }
1513
1514 if (nt != 0) {
1515 if (only_divides) { r=poly(BHEAD 1); ri=r[0]; break; }
1516 for (int i=0; i<AN.poly_num_vars; i++)
1517 r[ri+1+i] = 0;
1518 r[ri+1+var] = pow;
1519
1520 for (int i=0; i<ABS(nt); i++)
1521 r[ri+1+AN.poly_num_vars+i] = t[i];
1522
1523 r[ri] = 2+AN.poly_num_vars+ABS(nt);
1524 ri += r[ri];
1525 r[ri-1] = nt;
1526 }
1527 }
1528 }
1529
1530 q[0] = qi;
1531 r[0] = ri;
1532
1533 NumberFree(s,"poly::div_univar");
1534 NumberFree(t,"poly::div_univar");
1535
1536 if (q.modp!=0 || ltbinv!=NULL) NumberFree(ltbinv,"poly::div_univar");
1537 }
1538
1539 /*
1540 #] divmod_univar :
1541 #[ divmod_heap :
1542 */
1543
1544 /** Division of polynomials with a heap
1545 *
1546 * Description
1547 * ===========
1548 * Divides two multivariate polynomials. The next element of the
1549 * quotient/remainder is efficiently determined by using a heap. If
1550 * the product of the maximum power in all variables is small, a
1551 * hash table is used to add equal terms for extra speed.
1552 *
1553 * If the input flag check_div is set then if the result of any
1554 * coefficient division results in a non-zero remainder (indicating
1555 * that division has failed over the integers) the output flag div_fail
1556 * will be set and the division will be terminated early (q, r
1557 * will be incorrect). If check_div is not set then non-zero remainders
1558 * from coefficient division will be written into r.
1559 *
1560 * A heap element h is formatted as follows:
1561 * - h[0] = index in a
1562 * - h[1] = index in b
1563 * - h[2] = -1 (no hash is used)
1564 * - h[3] = length of coefficient with sign
1565 * - h[4...4+AN.poly_num_vars-1] = powers
1566 * - h[4+AN.poly_num_vars...4+h[3]-1] = coefficient
1567 *
1568 * Note: the hashing trick as in multiplication cannot be used
1569 * easily, since there is no tight upperbound on the exponents in
1570 * the answer.
1571
1572 * For details, see M. Monagan, "Polynomial Division using Dynamic
1573 * Array, Heaps, and Packed Exponent Vectors"
1574 */
divmod_heap(const poly & a,const poly & b,poly & q,poly & r,bool only_divides,bool check_div,bool & div_fail)1575 void poly::divmod_heap (const poly &a, const poly &b, poly &q, poly &r, bool only_divides, bool check_div, bool &div_fail) {
1576
1577 POLY_GETIDENTITY(a);
1578
1579 div_fail = false;
1580 q[0] = r[0] = 1;
1581
1582 WORD nmodq=0;
1583 UWORD *modq=NULL;
1584
1585 WORD nltbinv=0;
1586 UWORD *ltbinv=NULL;
1587 LONG oldpWorkPointer = AT.pWorkPointer;
1588
1589 bool both_mod_small=false;
1590
1591 if (q.modp!=0) {
1592 if (q.modn == 1) {
1593 modq = (UWORD *)&q.modp;
1594 nmodq = 1;
1595 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
1596 both_mod_small=true;
1597 }
1598 else {
1599 RaisPowCached(BHEAD q.modp,q.modn,&modq,&nmodq);
1600 }
1601 ltbinv = NumberMalloc("poly::div_heap-a");
1602
1603 if (both_mod_small) {
1604 WORD ltb = b[b[1]]*b[2+AN.poly_num_vars];
1605 GetModInverses(ltb + (ltb<0?q.modp:0), q.modp, (WORD*)ltbinv, NULL);
1606 nltbinv = 1;
1607 }
1608 else
1609 GetLongModInverses(BHEAD (UWORD *)&b[2+AN.poly_num_vars], b[b[1]], modq, nmodq, ltbinv, &nltbinv, NULL, NULL);
1610 }
1611
1612 // allocate heap
1613 int nb=b.number_of_terms();
1614 WantAddPointers(nb);
1615 WORD **heap = AT.pWorkSpace + AT.pWorkPointer;
1616 AT.pWorkPointer += nb;
1617
1618 for (int i=0; i<nb; i++)
1619 heap[i] = (WORD *) NumberMalloc("poly::div_heap-b");
1620 int nheap = 1;
1621 heap[0][0] = 1;
1622 heap[0][1] = 0;
1623 heap[0][2] = -1;
1624 WCOPY(&heap[0][3], &a[1], a[1]);
1625 heap[0][3] = a[a[1]];
1626
1627 int qi=1, ri=1;
1628
1629 int s = nb;
1630 WORD *t = (WORD *) NumberMalloc("poly::div_heap-c");
1631
1632 // insert contains element that still have to be inserted to the heap
1633 // (exists to avoid code duplication).
1634 vector<pair<int,int> > insert;
1635
1636 while (insert.size()>0 || nheap>0) {
1637
1638 q.check_memory(qi);
1639 r.check_memory(ri);
1640
1641 // collect a term t for the quotient/remainder
1642 t[0] = -1;
1643
1644 do {
1645
1646 WORD *p = heap[nheap];
1647 bool this_insert;
1648
1649 if (insert.empty()) {
1650 // extract element from the heap and prepare adding new ones
1651 this_insert = false;
1652
1653 pop_heap(BHEAD heap, nheap--);
1654 p = heap[nheap];
1655
1656 if (t[0] == -1) {
1657 WCOPY(t, p, (5+ABS(p[3])+AN.poly_num_vars));
1658 }
1659 else {
1660 // if both polynomials are modulo p^1, use integer calculus
1661 if (both_mod_small) {
1662 t[4+AN.poly_num_vars] = ((LONG)t[4+AN.poly_num_vars]*t[3] + p[4+AN.poly_num_vars]*p[3]) % q.modp;
1663 if (t[4+AN.poly_num_vars]==0)
1664 t[3]=0;
1665 else {
1666 if (t[4+AN.poly_num_vars] > +q.modp/2) t[4+AN.poly_num_vars] -= q.modp;
1667 if (t[4+AN.poly_num_vars] < -q.modp/2) t[4+AN.poly_num_vars] += q.modp;
1668 t[3] = SGN(t[4+AN.poly_num_vars]);
1669 t[4+AN.poly_num_vars] = ABS(t[4+AN.poly_num_vars]);
1670 }
1671 }
1672 else {
1673 // otherwise, use form long calculus
1674 AddLong ((UWORD *)&p[4+AN.poly_num_vars], p[3],
1675 (UWORD *)&t[4+AN.poly_num_vars], t[3],
1676 (UWORD *)&t[4+AN.poly_num_vars], &t[3]);
1677 if (q.modp!=0) TakeNormalModulus((UWORD *)&t[4+AN.poly_num_vars], &t[3],
1678 modq, nmodq, NOUNPACK);
1679 }
1680 }
1681 }
1682 else {
1683 // prepare adding an element of insert to the heap
1684 this_insert = true;
1685
1686 p[0] = insert.back().first;
1687 p[1] = insert.back().second;
1688 insert.pop_back();
1689 }
1690
1691 // add elements to the heap
1692 while (true) {
1693 // prepare the element
1694 if (p[1]==0) {
1695 p[0] += a[p[0]];
1696 if (p[0]==a[0]) break;
1697 WCOPY(&p[3], &a[p[0]], a[p[0]]);
1698 p[3] = p[2+p[3]];
1699 }
1700 else {
1701 if (!this_insert)
1702 p[1] += q[p[1]];
1703 this_insert = false;
1704
1705 if (p[1]==qi) { s++; break; }
1706
1707 for (int i=0; i<AN.poly_num_vars; i++)
1708 p[4+i] = b[p[0]+1+i] + q[p[1]+1+i];
1709
1710 // if both polynomials are modulo p^1, use integer calculus
1711 if (both_mod_small) {
1712 p[4+AN.poly_num_vars] = ((LONG)b[p[0]+1+AN.poly_num_vars]*b[p[0]+b[p[0]]-1]*
1713 q[p[1]+1+AN.poly_num_vars]*q[p[1]+q[p[1]]-1]) % q.modp;
1714 if (p[4+AN.poly_num_vars]==0)
1715 p[3]=0;
1716 else {
1717 if (p[4+AN.poly_num_vars] > +q.modp/2) p[4+AN.poly_num_vars] -= q.modp;
1718 if (p[4+AN.poly_num_vars] < -q.modp/2) p[4+AN.poly_num_vars] += q.modp;
1719 p[3] = SGN(p[4+AN.poly_num_vars]);
1720 p[4+AN.poly_num_vars] = ABS(p[4+AN.poly_num_vars]);
1721 }
1722 }
1723 else {
1724 // otherwise, use form long calculus
1725 MulLong((UWORD *)&b[p[0]+1+AN.poly_num_vars], b[p[0]+b[p[0]]-1],
1726 (UWORD *)&q[p[1]+1+AN.poly_num_vars], q[p[1]+q[p[1]]-1],
1727 (UWORD *)&p[4+AN.poly_num_vars], &p[3]);
1728 if (q.modp!=0) TakeNormalModulus((UWORD *)&p[4+AN.poly_num_vars], &p[3],
1729 modq, nmodq, NOUNPACK);
1730 }
1731
1732 p[3] *= -1;
1733 }
1734
1735 // no hashing
1736 p[2] = -1;
1737
1738 // add it to a heap element
1739 swap (heap[nheap],p);
1740 push_heap(BHEAD heap, ++nheap);
1741 break;
1742 }
1743 }
1744 while (t[0]==-1 || (nheap>0 && monomial_compare(BHEAD heap[0]+3, t+3)==0));
1745
1746 if (t[3] == 0) continue;
1747
1748 // check divisibility
1749 bool div = true;
1750 for (int i=0; i<AN.poly_num_vars; i++)
1751 if (t[4+i] < b[2+i]) div=false;
1752
1753 if (!div) {
1754 // not divisible, so append it to the remainder
1755 if (only_divides) { r=poly(BHEAD 1); ri=r[0]; break; }
1756 t[4 + AN.poly_num_vars + ABS(t[3])] = t[3];
1757 t[3] = 2 + AN.poly_num_vars + ABS(t[3]);
1758 r.termscopy(&t[3], ri, t[3]);
1759 ri += t[3];
1760 }
1761 else {
1762 // divisable, so divide coefficient as well
1763 WORD nq, nr;
1764
1765 if (q.modp==0) {
1766 DivLong((UWORD *)&t[4+AN.poly_num_vars], t[3],
1767 (UWORD *)&b[2+AN.poly_num_vars], b[b[1]],
1768 (UWORD *)&q[qi+1+AN.poly_num_vars], &nq,
1769 (UWORD *)&r[ri+1+AN.poly_num_vars], &nr);
1770 if(check_div && nr != 0)
1771 {
1772 // non-zero remainder from coefficient division => result not expressible in integers
1773 div_fail = true;
1774 break;
1775 }
1776 }
1777 else {
1778 // if both polynomials are modulo p^1, use integer calculus
1779 if (both_mod_small) {
1780 q[qi+1+AN.poly_num_vars] = ((LONG)t[4+AN.poly_num_vars]*t[3]*ltbinv[0]*nltbinv) % q.modp;
1781 if (q[qi+1+AN.poly_num_vars]==0)
1782 nq=0;
1783 else {
1784 if (q[qi+1+AN.poly_num_vars] > +q.modp/2) q[qi+1+AN.poly_num_vars] -= q.modp;
1785 if (q[qi+1+AN.poly_num_vars] < -q.modp/2) q[qi+1+AN.poly_num_vars] += q.modp;
1786 nq = SGN(q[qi+1+AN.poly_num_vars]);
1787 q[qi+1+AN.poly_num_vars] = ABS(q[qi+1+AN.poly_num_vars]);
1788 }
1789 }
1790 else {
1791 // otherwise, use form long calculus
1792 MulLong((UWORD *)&t[4+AN.poly_num_vars], t[3], ltbinv, nltbinv, (UWORD *)&q[qi+1+AN.poly_num_vars], &nq);
1793 TakeNormalModulus((UWORD *)&q[qi+1+AN.poly_num_vars], &nq, modq, nmodq, NOUNPACK);
1794 }
1795
1796 nr=0;
1797 }
1798
1799 // add terms to quotient and remainder
1800 if (nq != 0) {
1801 int bi = 1;
1802 for (int j=1; j<s; j++) {
1803 bi += b[bi];
1804 insert.push_back(make_pair(bi,qi));
1805 }
1806 s=1;
1807
1808 q[qi] = 2+AN.poly_num_vars+ABS(nq);
1809 for (int i=0; i<AN.poly_num_vars; i++)
1810 q[qi+1+i] = t[4+i] - b[2+i];
1811 qi += q[qi];
1812 q[qi-1] = nq;
1813 }
1814
1815 if (nr != 0) {
1816 r[ri] = 2+AN.poly_num_vars+ABS(nr);
1817 for (int i=0; i<AN.poly_num_vars; i++)
1818 r[ri+1+i] = t[4+i];
1819 ri += r[ri];
1820 r[ri-1] = nr;
1821 }
1822 }
1823 }
1824
1825 q[0] = qi;
1826 r[0] = ri;
1827
1828 for (int i=0; i<nb; i++)
1829 NumberFree(heap[i],"poly::div_heap-b");
1830
1831 NumberFree(t,"poly::div_heap-c");
1832
1833 if (q.modp!=0||ltbinv!=NULL) NumberFree(ltbinv,"poly::div_heap-a");
1834 AT.pWorkPointer = oldpWorkPointer;
1835 }
1836
1837 /*
1838 #] divmod_heap :
1839 #[ divmod :
1840 */
1841
1842 /** Polynomial division
1843 *
1844 * Description
1845 * ===========
1846 * This routine determines which division routine to use for
1847 * dividing two polynomials. The logic is as follows:
1848 * - If b consists of only one term, call divmod_one_term;
1849 * - Otherwise, if both are univariate and dense, call divmod_univar;
1850 * - Otherwise, call divmod_heap.
1851 */
divmod(const poly & a,const poly & b,poly & q,poly & r,bool only_divides)1852 void poly::divmod (const poly &a, const poly &b, poly &q, poly &r, bool only_divides) {
1853
1854 q.modp = r.modp = a.modp;
1855 q.modn = r.modn = a.modn;
1856
1857 if (a.is_zero()) {
1858 q[0]=1;
1859 r[0]=1;
1860 return;
1861 }
1862 if (b.is_zero()) {
1863 MLOCK(ErrorMessageLock);
1864 MesPrint ((char*)"ERROR: division by zero polynomial");
1865 MUNLOCK(ErrorMessageLock);
1866 Terminate(1);
1867 }
1868 if (b.is_one()) {
1869 q.check_memory(a[0]);
1870 q.termscopy(a.terms,0,a[0]);
1871 r[0]=1;
1872 return;
1873 }
1874
1875 if (b.is_monomial()) {
1876 divmod_one_term(a,b,q,r,only_divides);
1877 return;
1878 }
1879
1880 int vara = a.is_dense_univariate();
1881 int varb = b.is_dense_univariate();
1882
1883 if (vara!=-2 && varb!=-2 && (vara==-1 || varb==-1 || vara==varb)) {
1884 divmod_univar(a,b,q,r,MaX(vara,varb),only_divides);
1885 }
1886 else {
1887 bool div_fail;
1888 divmod_heap(a,b,q,r,only_divides,false,div_fail);
1889 }
1890 }
1891
1892 /*
1893 #] divmod :
1894 #[ divides :
1895 */
1896
1897 // returns whether a exactly divides b
divides(const poly & a,const poly & b)1898 bool poly::divides (const poly &a, const poly &b) {
1899 POLY_GETIDENTITY(a);
1900 poly q(BHEAD 0), r(BHEAD 0);
1901 divmod(b,a,q,r,true);
1902 return r.is_zero();
1903 }
1904
1905 /*
1906 #] divides :
1907 #[ div :
1908 */
1909
1910 // the quotient of two polynomials
div(const poly & a,const poly & b,poly & c)1911 void poly::div (const poly &a, const poly &b, poly &c) {
1912 POLY_GETIDENTITY(a);
1913 poly d(BHEAD 0);
1914 divmod(a,b,c,d,false);
1915 }
1916
1917 /*
1918 #] div :
1919 #[ mod :
1920 */
1921
1922 // the remainder of division of two polynomials
mod(const poly & a,const poly & b,poly & c)1923 void poly::mod (const poly &a, const poly &b, poly &c) {
1924 POLY_GETIDENTITY(a);
1925 poly d(BHEAD 0);
1926 divmod(a,b,d,c,false);
1927 }
1928
1929 /*
1930 #] mod :
1931 #[ copy operator :
1932 */
1933
1934 // copy operator
operator =(const poly & a)1935 poly & poly::operator= (const poly &a) {
1936
1937 if (&a != this) {
1938 modp = a.modp;
1939 modn = a.modn;
1940 check_memory(a[0]);
1941 termscopy(a.terms,0,a[0]);
1942 }
1943
1944 return *this;
1945 }
1946
1947 /*
1948 #] copy operator :
1949 #[ operator overloads :
1950 */
1951
1952 // binary operators for polynomial arithmetic
operator +(const poly & a) const1953 const poly poly::operator+ (const poly &a) const { POLY_GETIDENTITY(a); poly b(BHEAD 0); add(*this,a,b); return b; }
operator -(const poly & a) const1954 const poly poly::operator- (const poly &a) const { POLY_GETIDENTITY(a); poly b(BHEAD 0); sub(*this,a,b); return b; }
operator *(const poly & a) const1955 const poly poly::operator* (const poly &a) const { POLY_GETIDENTITY(a); poly b(BHEAD 0); mul(*this,a,b); return b; }
operator /(const poly & a) const1956 const poly poly::operator/ (const poly &a) const { POLY_GETIDENTITY(a); poly b(BHEAD 0); div(*this,a,b); return b; }
operator %(const poly & a) const1957 const poly poly::operator% (const poly &a) const { POLY_GETIDENTITY(a); poly b(BHEAD 0); mod(*this,a,b); return b; }
1958
1959 // assignment operators for polynomial arithmetic
operator +=(const poly & a)1960 poly& poly::operator+= (const poly &a) { return *this = *this + a; }
operator -=(const poly & a)1961 poly& poly::operator-= (const poly &a) { return *this = *this - a; }
operator *=(const poly & a)1962 poly& poly::operator*= (const poly &a) { return *this = *this * a; }
operator /=(const poly & a)1963 poly& poly::operator/= (const poly &a) { return *this = *this / a; }
operator %=(const poly & a)1964 poly& poly::operator%= (const poly &a) { return *this = *this % a; }
1965
1966 // comparison operators
operator ==(const poly & a) const1967 bool poly::operator== (const poly &a) const {
1968 for (int i=0; i<terms[0]; i++)
1969 if (terms[i] != a[i]) return 0;
1970 return 1;
1971 }
1972
operator !=(const poly & a) const1973 bool poly::operator!= (const poly &a) const { return !(*this == a); }
1974
1975 /*
1976 #] operator overloads :
1977 #[ num_terms :
1978 */
1979
number_of_terms() const1980 int poly::number_of_terms () const {
1981
1982 int res=0;
1983 for (int i=1; i<terms[0]; i+=terms[i])
1984 res++;
1985 return res;
1986 }
1987
1988 /*
1989 #] num_terms :
1990 #[ first_variable :
1991 */
1992
1993 // returns the lexcicographically first variable of a polynomial
first_variable() const1994 int poly::first_variable () const {
1995
1996 POLY_GETIDENTITY(*this);
1997
1998 int var = AN.poly_num_vars;
1999 for (int j=0; j<var; j++)
2000 if (terms[2+j]>0) return j;
2001 return var;
2002 }
2003
2004 /*
2005 #] first_variable :
2006 #[ all_variables :
2007 */
2008
2009 // returns a list of all variables of a polynomial
all_variables() const2010 const vector<int> poly::all_variables () const {
2011
2012 POLY_GETIDENTITY(*this);
2013
2014 vector<bool> used(AN.poly_num_vars, false);
2015
2016 for (int i=1; i<terms[0]; i+=terms[i])
2017 for (int j=0; j<AN.poly_num_vars; j++)
2018 if (terms[i+1+j]>0) used[j] = true;
2019
2020 vector<int> vars;
2021 for (int i=0; i<AN.poly_num_vars; i++)
2022 if (used[i]) vars.push_back(i);
2023
2024 return vars;
2025 }
2026
2027 /*
2028 #] all_variables :
2029 #[ sign :
2030 */
2031
2032 // returns the sign of the leading coefficient
sign() const2033 int poly::sign () const {
2034 if (terms[0]==1) return 0;
2035 return terms[terms[1]] > 0 ? 1 : -1;
2036 }
2037
2038 /*
2039 #] sign :
2040 #[ degree :
2041 */
2042
2043 // returns the degree of x of a polynomial (deg=-1 iff a=0)
degree(int x) const2044 int poly::degree (int x) const {
2045 int deg = -1;
2046 for (int i=1; i<terms[0]; i+=terms[i])
2047 deg = MaX(deg, terms[i+1+x]);
2048 return deg;
2049 }
2050
2051 /*
2052 #] degree :
2053 #[ total_degree :
2054 */
2055
2056 // returns the total degree of a polynomial (deg=-1 iff a=0)
total_degree() const2057 int poly::total_degree () const {
2058
2059 POLY_GETIDENTITY(*this);
2060
2061 int tot_deg = -1;
2062 for (int i=1; i<terms[0]; i+=terms[i]) {
2063 int deg=0;
2064 for (int j=0; j<AN.poly_num_vars; j++)
2065 deg += terms[i+1+j];
2066 tot_deg = MaX(deg, tot_deg);
2067 }
2068 return tot_deg;
2069 }
2070
2071 /*
2072 #] total_degree :
2073 #[ integer_lcoeff :
2074 */
2075
2076 // returns the integer coefficient of the leading monomial
integer_lcoeff() const2077 const poly poly::integer_lcoeff () const {
2078
2079 POLY_GETIDENTITY(*this);
2080
2081 poly res(BHEAD 0);
2082 res.modp = modp;
2083 res.modn = modn;
2084
2085 res.termscopy(&terms[1],1,terms[1]);
2086 res[0] = res[1] + 1; // length
2087 for (int i=0; i<AN.poly_num_vars; i++)
2088 res[2+i] = 0; // powers
2089
2090 return res;
2091 }
2092
2093 /*
2094 #] integer_lcoeff :
2095 #[ coefficient :
2096 */
2097
2098 // returns the polynomial coefficient of x^n
coefficient(int x,int n) const2099 const poly poly::coefficient (int x, int n) const {
2100
2101 POLY_GETIDENTITY(*this);
2102
2103 poly res(BHEAD 0);
2104 res.modp = modp;
2105 res.modn = modn;
2106 res[0] = 1;
2107
2108 for (int i=1; i<terms[0]; i+=terms[i])
2109 if (terms[i+1+x] == n) {
2110 res.check_memory(res[0]+terms[i]);
2111 res.termscopy(&terms[i], res[0], terms[i]);
2112 res[res[0]+1+x] -= n; // power of x
2113 res[0] += res[res[0]]; // length
2114 }
2115
2116 return res;
2117 }
2118
2119 /*
2120 #] coefficient :
2121 #[ lcoeff_multivar :
2122 */
2123
2124 // returns the leading coefficient with the polynomial viewed as a
2125 // polynomial in x, so that the result is a polynomial in the
2126 // remaining variables
lcoeff_univar(int x) const2127 const poly poly::lcoeff_univar (int x) const {
2128 return coefficient(x, degree(x));
2129 }
2130
2131 // returns the leading coefficient with the polynomial viewed as a
2132 // polynomial with coefficients in ZZ[x], so that the result
2133 // polynomial is in ZZ[x] too
lcoeff_multivar(int x) const2134 const poly poly::lcoeff_multivar (int x) const {
2135
2136 POLY_GETIDENTITY(*this);
2137
2138 poly res(BHEAD 0,modp,modn);
2139
2140 for (int i=1; i<terms[0]; i+=terms[i]) {
2141 bool same_powers = true;
2142 for (int j=0; j<AN.poly_num_vars; j++)
2143 if (j!=x && terms[i+1+j]!=terms[2+j]) {
2144 same_powers = false;
2145 break;
2146 }
2147 if (!same_powers) break;
2148
2149 res.check_memory(res[0]+terms[i]);
2150 res.termscopy(&terms[i],res[0],terms[i]);
2151 for (int k=0; k<AN.poly_num_vars; k++)
2152 if (k!=x) res[res[0]+1+k]=0;
2153
2154 res[0] += terms[i];
2155 }
2156
2157 return res;
2158 }
2159
2160 /*
2161 #] lcoeff_multivar :
2162 #[ derivative :
2163 */
2164
2165 // returns the derivative with respect to x
derivative(int x) const2166 const poly poly::derivative (int x) const {
2167
2168 POLY_GETIDENTITY(*this);
2169
2170 poly b(BHEAD 0);
2171 int bi=1;
2172
2173 for (int i=1; i<terms[0]; i+=terms[i]) {
2174
2175 int power = terms[i+1+x];
2176
2177 if (power > 0) {
2178 b.check_memory(bi);
2179 b.termscopy(&terms[i], bi, terms[i]);
2180 b[bi+1+x]--;
2181
2182 WORD nb = b[bi+b[bi]-1];
2183 Product((UWORD *)&b[bi+1+AN.poly_num_vars], &nb, power);
2184
2185 b[bi] = 2 + AN.poly_num_vars + ABS(nb);
2186 b[bi+b[bi]-1] = nb;
2187
2188 bi += b[bi];
2189 }
2190 }
2191
2192 b[0] = bi;
2193 b.setmod(modp, modn);
2194 return b;
2195 }
2196
2197 /*
2198 #] derivative :
2199 #[ is_zero :
2200 */
2201
2202 // returns whether the polynomial is zero
is_zero() const2203 bool poly::is_zero () const {
2204 return terms[0] == 1;
2205 }
2206
2207 /*
2208 #] is_zero :
2209 #[ is_one :
2210 */
2211
2212 // returns whether the polynomial is one
is_one() const2213 bool poly::is_one () const {
2214
2215 POLY_GETIDENTITY(*this);
2216
2217 if (terms[0] != 4+AN.poly_num_vars) return false;
2218 if (terms[1] != 3+AN.poly_num_vars) return false;
2219 for (int i=0; i<AN.poly_num_vars; i++)
2220 if (terms[2+i] != 0) return false;
2221 if (terms[2+AN.poly_num_vars]!=1) return false;
2222 if (terms[3+AN.poly_num_vars]!=1) return false;
2223
2224 return true;
2225 }
2226
2227 /*
2228 #] is_one :
2229 #[ is_integer :
2230 */
2231
2232 // returns whether the polynomial is an integer
is_integer() const2233 bool poly::is_integer () const {
2234
2235 POLY_GETIDENTITY(*this);
2236
2237 if (terms[0] == 1) return true;
2238 if (terms[0] != terms[1]+1) return false;
2239
2240 for (int j=0; j<AN.poly_num_vars; j++)
2241 if (terms[2+j] != 0)
2242 return false;
2243
2244 return true;
2245 }
2246
2247 /*
2248 #] is_integer :
2249 #[ is_monomial :
2250 */
2251
2252 // returns whether the polynomial consist of one term
is_monomial() const2253 bool poly::is_monomial () const {
2254 return terms[0]>1 && terms[0]==terms[1]+1;
2255 }
2256
2257 /*
2258 #] is_monomial :
2259 #[ is_dense_univariate :
2260 */
2261
2262 /** Dense univariate detection
2263 *
2264 * Description
2265 * ===========
2266 * This method returns whether the polynomial is dense and
2267 * univariate. The possible return values are:
2268 *
2269 * -2 is not dense univariate
2270 * -1 is no variables
2271 * n>=0 is univariate in n
2272 *
2273 * Notes
2274 * =====
2275 * A univariate polynomial is considered dense iff more than half of
2276 * the coefficients a_0...a_deg are non-zero.
2277 */
is_dense_univariate() const2278 int poly::is_dense_univariate () const {
2279
2280 POLY_GETIDENTITY(*this);
2281
2282 int num_terms=0, res=-1;
2283
2284 // test univariate
2285 for (int i=1; i<terms[0]; i+=terms[i]) {
2286 for (int j=0; j<AN.poly_num_vars; j++)
2287 if (terms[i+1+j] > 0) {
2288 if (res == -1) res = j;
2289 if (res != j) return -2;
2290 }
2291
2292 num_terms++;
2293 }
2294
2295 // constant polynomial
2296 if (res == -1) return -1;
2297
2298 // test density
2299 int deg = terms[2+res];
2300 if (2*num_terms < deg+1) return -2;
2301
2302 return res;
2303 }
2304
2305 /*
2306 #] is_dense_univariate :
2307 #[ simple_poly (small) :
2308 */
2309
2310 // returns the polynomial (x-a)^b mod p^n with a small
simple_poly(PHEAD int x,int a,int b,int p,int n)2311 const poly poly::simple_poly (PHEAD int x, int a, int b, int p, int n) {
2312
2313 poly tmp(BHEAD 0,p,n);
2314
2315 int idx=1;
2316 tmp[idx++] = 3 + AN.poly_num_vars; // length
2317 for (int i=0; i<AN.poly_num_vars; i++)
2318 tmp[idx++] = i==x ? 1 : 0; // powers
2319 tmp[idx++] = 1; // coefficient
2320 tmp[idx++] = 1; // length coefficient
2321
2322 if (a != 0) {
2323 tmp[idx++] = 3 + AN.poly_num_vars; // length
2324 for (int i=0; i<AN.poly_num_vars; i++) tmp[idx++] = 0; // powers
2325 tmp[idx++] = ABS(a); // coefficient
2326 tmp[idx++] = -SGN(a); // length coefficient
2327 }
2328
2329 tmp[0] = idx; // length
2330
2331 if (b == 1) return tmp;
2332
2333 poly res(BHEAD 1,p,n);
2334
2335 while (b!=0) {
2336 if (b&1) res*=tmp;
2337 tmp*=tmp;
2338 b>>=1;
2339 }
2340
2341 return res;
2342 }
2343
2344 /*
2345 #] simple_poly (small) :
2346 #[ simple_poly (large) :
2347 */
2348
2349 // returns the polynomial (x-a)^b mod p^n with a large
simple_poly(PHEAD int x,const poly & a,int b,int p,int n)2350 const poly poly::simple_poly (PHEAD int x, const poly &a, int b, int p, int n) {
2351
2352 poly res(BHEAD 1,p,n);
2353 poly tmp(BHEAD 0,p,n);
2354
2355 int idx=1;
2356
2357 tmp[idx++] = 3 + AN.poly_num_vars; // length
2358 for (int i=0; i<AN.poly_num_vars; i++)
2359 tmp[idx++] = i==x ? 1 : 0; // powers
2360 tmp[idx++] = 1; // coefficient
2361 tmp[idx++] = 1; // length coefficient
2362
2363 if (!a.is_zero()) {
2364 tmp[idx++] = 2 + AN.poly_num_vars + ABS(a[a[0]-1]); // length
2365 for (int i=0; i<AN.poly_num_vars; i++) tmp[idx++] = 0; // powers
2366 for (int i=0; i<ABS(a[a[0]-1]); i++)
2367 tmp[idx++] = a[2 + AN.poly_num_vars + i]; // coefficient
2368 tmp[idx++] = -a[a[0]-1]; // length coefficient
2369 }
2370
2371 tmp[0] = idx; // length
2372
2373 while (b!=0) {
2374 if (b&1) res*=tmp;
2375 tmp*=tmp;
2376 b>>=1;
2377 }
2378
2379 return res;
2380 }
2381
2382 /*
2383 #] simple_poly (large) :
2384 #[ get_variables :
2385 */
2386
2387 // gets all variables in the expressions and stores them in AN.poly_vars
2388 // (it is assumed that AN.poly_vars=NULL)
get_variables(PHEAD vector<WORD * > es,bool with_arghead,bool sort_vars)2389 void poly::get_variables (PHEAD vector<WORD *> es, bool with_arghead, bool sort_vars) {
2390
2391 AN.poly_num_vars = 0;
2392
2393 vector<int> vars;
2394 vector<int> degrees;
2395 map<int,int> var_to_idx;
2396
2397 // extract all variables
2398 for (int ei=0; ei<(int)es.size(); ei++) {
2399 WORD *e = es[ei];
2400
2401 // fast notation
2402 if (*e == -SNUMBER) {
2403 }
2404 else if (*e == -SYMBOL) {
2405 if (!var_to_idx.count(e[1])) {
2406 vars.push_back(e[1]);
2407 var_to_idx[e[1]] = AN.poly_num_vars++;
2408 degrees.push_back(1);
2409 }
2410 }
2411 else {
2412 for (int i=with_arghead ? ARGHEAD : 0; with_arghead ? i<e[0] : e[i]!=0; i+=e[i]) {
2413 if (i+1<i+e[i]-ABS(e[i+e[i]-1]) && e[i+1]!=SYMBOL) {
2414 MLOCK(ErrorMessageLock);
2415 MesPrint ((char*)"ERROR: polynomials and polyratfuns must contain symbols only");
2416 MUNLOCK(ErrorMessageLock);
2417 Terminate(1);
2418 }
2419
2420 for (int j=i+3; j<i+e[i]-ABS(e[i+e[i]-1]); j+=2)
2421 if (!var_to_idx.count(e[j])) {
2422 vars.push_back(e[j]);
2423 var_to_idx[e[j]] = AN.poly_num_vars++;
2424 degrees.push_back(e[j+1]);
2425 }
2426 else {
2427 degrees[var_to_idx[e[j]]] = MaX(degrees[var_to_idx[e[j]]], e[j+1]);
2428 }
2429 }
2430 }
2431 }
2432
2433 // make sure an eventual FACTORSYMBOL appear as last
2434 if (var_to_idx.count(FACTORSYMBOL)) {
2435 int i = var_to_idx[FACTORSYMBOL];
2436 while (i+1<AN.poly_num_vars) {
2437 swap(vars[i], vars[i+1]);
2438 swap(degrees[i], degrees[i+1]);
2439 i++;
2440 }
2441 degrees[i] = -1; // makes sure it stays last
2442 }
2443
2444 // AN.poly_vars will be deleted in calling functions from polywrap.c
2445 // For that we use the function poly_free_poly_vars
2446 // We add one to the allocation to avoid copying in poly_divmod
2447
2448 if ( (AN.poly_num_vars+1)*sizeof(WORD) > (size_t)(AM.MaxTer) ) {
2449 // This can happen only in expressions with excessively many variables.
2450 AN.poly_vars = (WORD *)Malloc1((AN.poly_num_vars+1)*sizeof(WORD), "AN.poly_vars");
2451 AN.poly_vars_type = 1;
2452 }
2453 else {
2454 AN.poly_vars = TermMalloc("AN.poly_vars");
2455 AN.poly_vars_type = 0;
2456 }
2457
2458 for (int i=0; i<AN.poly_num_vars; i++)
2459 AN.poly_vars[i] = vars[i];
2460
2461 if (sort_vars) {
2462 // bubble sort variables in decreasing order of degree
2463 // (this seems better for factorization)
2464 for (int i=0; i<AN.poly_num_vars; i++)
2465 for (int j=0; j+1<AN.poly_num_vars; j++)
2466 if (degrees[j] < degrees[j+1]) {
2467 swap(degrees[j], degrees[j+1]);
2468 swap(AN.poly_vars[j], AN.poly_vars[j+1]);
2469 }
2470 }
2471 else {
2472 // sort lexicographically
2473 sort(AN.poly_vars, AN.poly_vars+AN.poly_num_vars - var_to_idx.count(FACTORSYMBOL));
2474 }
2475 }
2476
2477 /*
2478 #] get_variables :
2479 #[ argument_to_poly :
2480 */
2481
2482 // converts a form expression to a polynomial class "poly"
argument_to_poly(PHEAD WORD * e,bool with_arghead,bool sort_univar,poly * denpoly)2483 const poly poly::argument_to_poly (PHEAD WORD *e, bool with_arghead, bool sort_univar, poly *denpoly) {
2484
2485 map<int,int> var_to_idx;
2486 for (int i=0; i<AN.poly_num_vars; i++)
2487 var_to_idx[AN.poly_vars[i]] = i;
2488
2489 poly res(BHEAD 0);
2490
2491 // fast notation
2492 if (*e == -SNUMBER) {
2493 if (denpoly!=NULL) *denpoly = poly(BHEAD 1);
2494
2495 if (e[1] == 0) {
2496 res[0] = 1;
2497 return res;
2498 }
2499 else {
2500 res[0] = 4 + AN.poly_num_vars;
2501 res[1] = 3 + AN.poly_num_vars;
2502 for (int i=0; i<AN.poly_num_vars; i++)
2503 res[2+i] = 0;
2504 res[2+AN.poly_num_vars] = ABS(e[1]);
2505 res[3+AN.poly_num_vars] = SGN(e[1]);
2506
2507 return res;
2508 }
2509 }
2510
2511 if (*e == -SYMBOL) {
2512 if (denpoly!=NULL) *denpoly = poly(BHEAD 1);
2513
2514 res[0] = 4 + AN.poly_num_vars;
2515 res[1] = 3 + AN.poly_num_vars;
2516 for (int i=0; i<AN.poly_num_vars; i++)
2517 res[2+i] = 0;
2518 res[2+var_to_idx.find(e[1])->second] = 1;
2519 res[2+AN.poly_num_vars] = 1;
2520 res[3+AN.poly_num_vars] = 1;
2521 return res;
2522 }
2523
2524 // find LCM of denominators of all terms
2525 WORD nden=1, npro=0, ngcd=0, ndum=0;
2526 UWORD *den = NumberMalloc("poly::argument_to_poly");
2527 UWORD *pro = NumberMalloc("poly::argument_to_poly");
2528 UWORD *gcd = NumberMalloc("poly::argument_to_poly");
2529 UWORD *dum = NumberMalloc("poly::argument_to_poly");
2530 den[0]=1;
2531
2532 for (int i=with_arghead ? ARGHEAD : 0; with_arghead ? i<e[0] : e[i]!=0; i+=e[i]) {
2533 int ncoe = ABS(e[i+e[i]-1]/2);
2534 UWORD *coe = (UWORD *)&e[i+e[i]-ncoe-1];
2535 while (ncoe>0 && coe[ncoe-1]==0) ncoe--;
2536 MulLong(den,nden,coe,ncoe,pro,&npro);
2537 GcdLong(BHEAD den,nden,coe,ncoe,gcd,&ngcd);
2538 DivLong(pro,npro,gcd,ngcd,den,&nden,dum,&ndum);
2539 }
2540
2541 if (denpoly!=NULL) *denpoly = poly(BHEAD den, nden);
2542
2543 int ri=1;
2544
2545 // ordinary notation
2546 for (int i=with_arghead ? ARGHEAD : 0; with_arghead ? i<e[0] : e[i]!=0; i+=e[i]) {
2547 res.check_memory(ri);
2548 WORD nc = e[i+e[i]-1]; // length coefficient
2549 for (int j=0; j<AN.poly_num_vars; j++)
2550 res[ri+1+j]=0; // powers=0
2551 WCOPY(dum, &e[i+e[i]-ABS(nc)], ABS(nc)); // coefficient to dummy
2552 nc /= 2; // remove denominator
2553 Mully(BHEAD dum, &nc, den, nden); // multiply with overall den
2554 res.termscopy((WORD *)dum, ri+1+AN.poly_num_vars, ABS(nc)); // coefficient to res
2555 res[ri] = ABS(nc) + AN.poly_num_vars + 2; // length
2556 res[ri+res[ri]-1] = nc; // length coefficient
2557 for (int j=i+3; j<i+e[i]-ABS(e[i+e[i]-1]); j+=2)
2558 res[ri+1+var_to_idx.find(e[j])->second] = e[j+1]; // powers
2559 ri += res[ri]; // length
2560 }
2561
2562 res[0] = ri;
2563
2564 // normalize, since the Form order is probably not the polynomial order
2565 // for multiple variables
2566
2567 if (sort_univar || AN.poly_num_vars>1)
2568 res.normalize();
2569
2570 NumberFree(den,"poly::argument_to_poly");
2571 NumberFree(pro,"poly::argument_to_poly");
2572 NumberFree(gcd,"poly::argument_to_poly");
2573 NumberFree(dum,"poly::argument_to_poly");
2574
2575 return res;
2576 }
2577
2578 /*
2579 #] argument_to_poly :
2580 #[ poly_to_argument :
2581 */
2582
2583 // converts a polynomial class "poly" to a form expression
poly_to_argument(const poly & a,WORD * res,bool with_arghead)2584 void poly::poly_to_argument (const poly &a, WORD *res, bool with_arghead) {
2585
2586 POLY_GETIDENTITY(a);
2587
2588 // special case: a=0
2589 if (a[0]==1) {
2590 if (with_arghead) {
2591 res[0] = -SNUMBER;
2592 res[1] = 0;
2593 }
2594 else {
2595 res[0] = 0;
2596 }
2597 return;
2598 }
2599
2600 if (with_arghead) {
2601 res[1] = AN.poly_num_vars>1 ? DIRTYFLAG : 0; // dirty flag
2602 for (int i=2; i<ARGHEAD; i++)
2603 res[i] = 0; // remainder of arghead
2604 }
2605
2606 int L = with_arghead ? ARGHEAD : 0;
2607
2608 for (int i=1; i!=a[0]; i+=a[i]) {
2609
2610 res[L]=1; // length
2611
2612 bool first=true;
2613
2614 for (int j=0; j<AN.poly_num_vars; j++)
2615 if (a[i+1+j] > 0) {
2616 if (first) {
2617 first=false;
2618 res[L+1] = 1; // symbols
2619 res[L+2] = 2; // length
2620 }
2621 res[L+1+res[L+2]++] = AN.poly_vars[j]; // symbol
2622 res[L+1+res[L+2]++] = a[i+1+j]; // power
2623 }
2624
2625 if (!first) res[L] += res[L+2]; // fix length
2626
2627 WORD nc = a[i+a[i]-1];
2628 WCOPY(&res[L+res[L]], &a[i+a[i]-1-ABS(nc)], ABS(nc)); // numerator
2629
2630 res[L] += ABS(nc); // fix length
2631 memset(&res[L+res[L]], 0, ABS(nc)*sizeof(WORD)); // denominator one
2632 res[L+res[L]] = 1; // denominator one
2633 res[L] += ABS(nc); // fix length
2634 res[L+res[L]] = SGN(nc) * (2*ABS(nc)+1); // length of coefficient
2635 res[L]++; // fix length
2636 L += res[L]; // fix length
2637 }
2638
2639 if (with_arghead) {
2640 res[0] = L;
2641 // convert to fast notation if possible
2642 ToFast(res,res);
2643 }
2644 else {
2645 res[L] = 0;
2646 }
2647 }
2648
2649 /*
2650 #] poly_to_argument :
2651 #[ poly_to_argument_with_den :
2652 */
2653
2654 // converts a polynomial class "poly" divided by a number (nden, den) to a form expression
2655 // cf. poly::poly_to_argument()
poly_to_argument_with_den(const poly & a,WORD nden,const UWORD * den,WORD * res,bool with_arghead)2656 void poly::poly_to_argument_with_den (const poly &a, WORD nden, const UWORD *den, WORD *res, bool with_arghead) {
2657
2658 POLY_GETIDENTITY(a);
2659
2660 // special case: a=0
2661 if (a[0]==1) {
2662 if (with_arghead) {
2663 res[0] = -SNUMBER;
2664 res[1] = 0;
2665 }
2666 else {
2667 res[0] = 0;
2668 }
2669 return;
2670 }
2671
2672 if (with_arghead) {
2673 res[1] = AN.poly_num_vars>1 ? DIRTYFLAG : 0; // dirty flag
2674 for (int i=2; i<ARGHEAD; i++)
2675 res[i] = 0; // remainder of arghead
2676 }
2677
2678 WORD nden1;
2679 UWORD *den1 = (UWORD *)NumberMalloc("poly_to_argument_with_den");
2680
2681 int L = with_arghead ? ARGHEAD : 0;
2682
2683 for (int i=1; i!=a[0]; i+=a[i]) {
2684
2685 res[L]=1; // length
2686
2687 bool first=true;
2688
2689 for (int j=0; j<AN.poly_num_vars; j++)
2690 if (a[i+1+j] > 0) {
2691 if (first) {
2692 first=false;
2693 res[L+1] = 1; // symbols
2694 res[L+2] = 2; // length
2695 }
2696 res[L+1+res[L+2]++] = AN.poly_vars[j]; // symbol
2697 res[L+1+res[L+2]++] = a[i+1+j]; // power
2698 }
2699
2700 if (!first) res[L] += res[L+2]; // fix length
2701
2702 // numerator
2703 WORD nc = a[i+a[i]-1];
2704 WCOPY(&res[L+res[L]], &a[i+a[i]-1-ABS(nc)], ABS(nc));
2705
2706 // denominator
2707 nden1 = nden;
2708 WCOPY(den1, den, ABS(nden));
2709
2710 if (nden != 1 || den[0] != 1) {
2711 // remove gcd(num,den)
2712 Simplify(BHEAD (UWORD *)&res[L+res[L]], &nc, den1, &nden1);
2713 }
2714
2715 Pack((UWORD *)&res[L+res[L]], &nc, den1, nden1); // format
2716 res[L] += 2*ABS(nc)+1; // fix length
2717 res[L+res[L]-1] = SGN(nc)*(2*ABS(nc)+1); // length of coefficient
2718 L += res[L]; // fix length
2719 }
2720
2721 NumberFree(den1, "poly_to_argument_with_den");
2722
2723 if (with_arghead) {
2724 res[0] = L;
2725 // convert to fast notation if possible
2726 ToFast(res,res);
2727 }
2728 else {
2729 res[L] = 0;
2730 }
2731 }
2732
2733 /*
2734 #] poly_to_argument_with_den :
2735 #[ size_of_form_notation :
2736 */
2737
2738 // the size of the polynomial in form notation (without argheads and fast notation)
size_of_form_notation() const2739 int poly::size_of_form_notation () const {
2740
2741 POLY_GETIDENTITY(*this);
2742
2743 // special case: a=0
2744 if (terms[0]==1) return 0;
2745
2746 int len = 0;
2747
2748 for (int i=1; i!=terms[0]; i+=terms[i]) {
2749 len++;
2750 int npow = 0;
2751 for (int j=0; j<AN.poly_num_vars; j++)
2752 if (terms[i+1+j] > 0) npow++;
2753 if (npow > 0) len += 2*npow + 2;
2754 len += 2 * ABS(terms[i+terms[i]-1]) + 1;
2755 }
2756
2757 return len;
2758 }
2759
2760 /*
2761 #] size_of_form_notation :
2762 #[ size_of_form_notation_with_den :
2763 */
2764
2765 // the size of the polynomial divided by a number (its size is given by nden)
2766 // in form notation (without argheads and fast notation)
2767 // cf. poly::size_of_form_notation()
size_of_form_notation_with_den(WORD nden) const2768 int poly::size_of_form_notation_with_den (WORD nden) const {
2769
2770 POLY_GETIDENTITY(*this);
2771
2772 // special case: a=0
2773 if (terms[0]==1) return 0;
2774
2775 nden = ABS(nden);
2776 int len = 0;
2777
2778 for (int i=1; i!=terms[0]; i+=terms[i]) {
2779 len++;
2780 int npow = 0;
2781 for (int j=0; j<AN.poly_num_vars; j++)
2782 if (terms[i+1+j] > 0) npow++;
2783 if (npow > 0) len += 2*npow + 2;
2784 len += 2 * MaX(ABS(terms[i+terms[i]-1]), nden) + 1;
2785 }
2786
2787 return len;
2788 }
2789
2790 /*
2791 #] size_of_form_notation_with_den :
2792 #[ to_coefficient_list :
2793 */
2794
2795 // returns the coefficient list of a univariate polynomial
to_coefficient_list(const poly & a)2796 const vector<WORD> poly::to_coefficient_list (const poly &a) {
2797
2798 POLY_GETIDENTITY(a);
2799
2800 if (a.is_zero()) return vector<WORD>();
2801
2802 int x = a.first_variable();
2803 if (x == AN.poly_num_vars) x=0;
2804
2805 vector<WORD> res(1+a[2+x],0);
2806
2807 for (int i=1; i<a[0]; i+=a[i])
2808 res[a[i+1+x]] = (a[i+a[i]-1] * a[i+a[i]-2]) % a.modp;
2809
2810 return res;
2811 }
2812
2813 /*
2814 #] to_coefficient_list :
2815 #[ coefficient_list_divmod :
2816 */
2817
2818 // divides two polynomials represented by coefficient lists
coefficient_list_divmod(const vector<WORD> & a,const vector<WORD> & b,WORD p,int divmod)2819 const vector<WORD> poly::coefficient_list_divmod (const vector<WORD> &a, const vector<WORD> &b, WORD p, int divmod) {
2820
2821 int bsize = (int)b.size();
2822 while (b[bsize-1]==0) bsize--;
2823
2824 WORD inv;
2825 GetModInverses(b[bsize-1] + (b[bsize-1]<0?p:0), p, &inv, NULL);
2826
2827 vector<WORD> q(a.size(),0);
2828 vector<WORD> r(a);
2829
2830 while ((int)r.size() >= bsize) {
2831 LONG mul = ((LONG)inv * r.back()) % p;
2832 int off = r.size()-bsize;
2833 q[off] = mul;
2834 for (int i=0; i<bsize; i++)
2835 r[off+i] = (r[off+i] - mul*b[i]) % p;
2836 while (r.size()>0 && r.back()==0)
2837 r.pop_back();
2838 }
2839
2840 if (divmod==0) {
2841 while (q.size()>0 && q.back()==0)
2842 q.pop_back();
2843 return q;
2844 }
2845 else {
2846 while (r.size()>0 && r.back()==0)
2847 r.pop_back();
2848 return r;
2849 }
2850 }
2851
2852 /*
2853 #] coefficient_list_divmod :
2854 #[ from_coefficient_list :
2855 */
2856
2857 // converts a coefficient list to a "poly"
from_coefficient_list(PHEAD const vector<WORD> & a,int x,WORD p)2858 const poly poly::from_coefficient_list (PHEAD const vector<WORD> &a, int x, WORD p) {
2859
2860 poly res(BHEAD 0);
2861 int ri=1;
2862
2863 for (int i=(int)a.size()-1; i>=0; i--)
2864 if (a[i] != 0) {
2865 res.check_memory(ri);
2866 res[ri] = AN.poly_num_vars+3; // length
2867 for (int j=0; j<AN.poly_num_vars; j++)
2868 res[ri+1+j]=0; // powers
2869 res[ri+1+x] = i; // power of x
2870 res[ri+1+AN.poly_num_vars] = ABS(a[i]); // coefficient
2871 res[ri+1+AN.poly_num_vars+1] = SGN(a[i]); // sign
2872 ri += res[ri];
2873 }
2874
2875 res[0]=ri; // total length
2876 res.setmod(p,1);
2877
2878 return res;
2879 }
2880
2881 /*
2882 #] from_coefficient_list :
2883 */
2884