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