53poly::poly (PHEAD
int a, WORD modp, WORD modn):
54 size_of_terms(AM.MaxTer/(LONG)sizeof(WORD)),
60 terms = TermMalloc(
"poly::poly(int)");
66 terms[0] = 4 + AN.poly_num_vars;
67 terms[1] = 3 + AN.poly_num_vars;
68 for (
int i=0; i<AN.poly_num_vars; i++) terms[2+i] = 0;
69 terms[2+AN.poly_num_vars] = ABS(a);
70 terms[3+AN.poly_num_vars] = a>0 ? 1 : -1;
73 if (modp!=-1) setmod(modp,modn);
82poly::poly (PHEAD
const UWORD *a, WORD na, WORD modp, WORD modn):
83 size_of_terms(AM.MaxTer/(LONG)sizeof(WORD)),
89 terms = TermMalloc(
"poly::poly(UWORD*,WORD)");
92 while (*(a+ABS(na)-1)==0)
95 terms[0] = 3 + AN.poly_num_vars + ABS(na);
96 terms[1] = terms[0] - 1;
97 for (
int i=0; i<AN.poly_num_vars; i++) terms[2+i] = 0;
98 termscopy((WORD *)a, 2+AN.poly_num_vars, ABS(na));
99 terms[2+AN.poly_num_vars+ABS(na)] = na;
101 if (modp!=-1) setmod(modp,modn);
110poly::poly (
const poly &a, WORD modp, WORD modn):
111 size_of_terms(AM.MaxTer/(LONG)sizeof(WORD))
116 terms = TermMalloc(
"poly::poly(poly&)");
120 if (modp!=-1) setmod(modp,modn);
131 POLY_GETIDENTITY(*
this);
133 if (size_of_terms == AM.MaxTer/(LONG)
sizeof(WORD))
134 TermFree(terms,
"poly::~poly");
136 M_free(terms,
"poly::~poly");
145void poly::expand_memory (
int i) {
147 POLY_GETIDENTITY(*
this);
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);
157 WORD *newterms = (WORD *)Malloc1(new_size_of_terms *
sizeof(WORD),
"poly::expand_memory");
158 WCOPY(newterms, terms, size_of_terms);
160 if (size_of_terms == AM.MaxTer/(LONG)
sizeof(WORD))
161 TermFree(terms,
"poly::expand_memory");
163 M_free(terms,
"poly::expand_memory");
166 size_of_terms = new_size_of_terms;
175void poly::setmod(WORD _modp, WORD _modn) {
177 POLY_GETIDENTITY(*
this);
179 if (_modp>0 && (_modp!=modp || _modn<modn)) {
188 modq = (UWORD *)&modp;
197 coefficients_modulo(modq,nmodq,smallq);
211void poly::coefficients_modulo (UWORD *a, WORD na,
bool small) {
213 POLY_GETIDENTITY(*
this);
216 for (
int i=1, di; i<terms[0]; i+=di) {
220 for (
int k=0; k<di; k++)
221 terms[j+k] = terms[i+k];
223 WORD n = terms[j+terms[j]-1];
225 if (ABS(n)==1 && small) {
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)
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]);
238 TakeNormalModulus((UWORD *)&terms[j+1+AN.poly_num_vars], &n, a, na, NOUNPACK);
241 terms[j] = 2+AN.poly_num_vars+ABS(n);
242 terms[j+terms[j]-1] = n;
256const string int_to_string (WORD x) {
258 snprintf (res, 41,
"%i",x);
263const string poly::to_string()
const {
265 POLY_GETIDENTITY(*
this);
270 UBYTE *scoeff = (UBYTE *)NumberMalloc(
"poly::to_string");
276 for (
int i=1; i<terms[0]; i+=terms[i]) {
279 WORD ncoeff = terms[i+terms[i]-1];
288 if (ncoeff==1 && terms[i+terms[i]-1-ncoeff]==1) {
294 PrtLong((UWORD*)&terms[i+terms[i]-1-ncoeff], ncoeff, scoeff);
295 res += string((
char *)scoeff);
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]);
310 if (!printtimes) res +=
"1";
317 res += int_to_string(modp);
320 res += int_to_string(modn);
325 NumberFree(scoeff,
"poly::to_string");
336ostream& operator<< (ostream &out,
const poly &a) {
337 return out << a.to_string();
346int 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];
358const poly & poly::normalize() {
360 POLY_GETIDENTITY(*
this);
367 modq = (UWORD *)&modp;
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]))
386 for (
int i=1; i<terms[0]; i+=terms[i])
387 p[nterms++] = terms + i;
392 if (size_of_terms == AM.MaxTer/(LONG)
sizeof(WORD)) {
393 tmp = (WORD *)TermMalloc(
"poly::normalize");
396 tmp = (WORD *)Malloc1(size_of_terms *
sizeof(WORD),
"poly::normalize");
403 for (
int i=0; i<nterms; i++) {
404 if (i>0 && monomial_compare(BHEAD &tmp[j], p[i])==0) {
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);
411 tmp[j+1+AN.poly_num_vars+ABS(ncoeff)] = ncoeff;
412 tmp[j] = 2+AN.poly_num_vars+ABS(ncoeff);
418 WCOPY(&tmp[j],p[i],p[i][0]);
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;
431 if (tmp[j+tmp[j]-1]==0) {
440 WCOPY(terms,tmp,tmp[0]);
444 if (size_of_terms == AM.MaxTer/(LONG)
sizeof(WORD))
445 TermFree(tmp,
"poly::normalize");
447 M_free(tmp,
"poly::normalize");
449 AT.pWorkPointer = poffset;
461WORD poly::last_monomial_index ()
const {
462 POLY_GETIDENTITY(*
this);
463 return terms[0] - ABS(terms[terms[0]-1]) - AN.poly_num_vars - 2;
468WORD * poly::last_monomial ()
const {
469 return &terms[last_monomial_index()];
477int poly::compare_degree_vector(
const poly & b)
const {
478 POLY_GETIDENTITY(*
this);
481 if (terms[0] == 1 && b[0] == 1)
return 0;
482 if (terms[0] == 1)
return -1;
483 if (b[0] == 1)
return 1;
485 for (
int i = 0; i < AN.poly_num_vars; i++) {
487 int d2 = b.degree(i);
488 if (d1 != d2)
return d1 - d2;
499std::vector<int> poly::degree_vector()
const {
500 POLY_GETIDENTITY(*
this);
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);
515int poly::compare_degree_vector(
const vector<int> & b)
const {
516 POLY_GETIDENTITY(*
this);
518 if (terms[0] == 1)
return -1;
520 for (
int i = 0; i < AN.poly_num_vars; i++) {
522 if (d1 != b[i])
return d1 - b[i];
544 bool both_mod_small=
false;
548 modq = (UWORD *)&c.modp;
550 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
560 while (ai<a[0] || bi<b[0]) {
564 int cmp = ai<a[0] && bi<b[0] ? monomial_compare(BHEAD &a[ai],&b[bi]) : 0;
566 if (bi==b[0] || cmp>0) {
568 c.termscopy(&a[ai],ci,a[ai]);
572 else if (ai==a[0] || cmp<0) {
574 c.termscopy(&b[bi],ci,b[bi]);
580 c.termscopy(&a[ai],ci,MaX(a[ai],b[bi]));
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]);
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);
599 c[ci] = 2+AN.poly_num_vars+ABS(nc);
628 bool both_mod_small=
false;
632 modq = (UWORD *)&c.modp;
634 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
644 while (ai<a[0] || bi<b[0]) {
648 int cmp = ai<a[0] && bi<b[0] ? monomial_compare(BHEAD &a[ai],&b[bi]) : 0;
650 if (bi==b[0] || cmp>0) {
652 c.termscopy(&a[ai],ci,a[ai]);
656 else if (ai==a[0] || cmp<0) {
658 c.termscopy(&b[bi],ci,b[bi]);
665 c.termscopy(&a[ai],ci,MaX(a[ai],b[bi]));
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]);
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],
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);
684 c[ci] = 2+AN.poly_num_vars+ABS(nc);
703void poly::pop_heap (PHEAD WORD **heap,
int n) {
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)) {
713 if (monomial_compare(BHEAD heap[2*i+1]+3, heap[2*i+2]+3)>0) {
714 swap(heap[i], heap[2*i+1]);
718 swap(heap[i], heap[2*i+2]);
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]);
735void poly::push_heap (PHEAD WORD **heap,
int n) {
739 while (i>0 && monomial_compare(BHEAD heap[i]+3, heap[(i-1)/2]+3) > 0) {
740 swap(heap[(i-1)/2], heap[i]);
751void poly::mul_one_term (
const poly &a,
const poly &b,
poly &c) {
760 bool both_mod_small=
false;
764 modq = (UWORD *)&c.modp;
766 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
774 for (
int ai=1; ai<a[0]; ai+=a[ai])
775 for (
int bi=1; bi<b[0]; bi+=b[bi]) {
779 for (
int i=0; i<AN.poly_num_vars; i++)
780 c[ci+1+i] = a[ai+1+i] + b[bi+1+i];
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);
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);
800 c[ci] = 2+AN.poly_num_vars+ABS(nc);
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]);
825void poly::mul_univar (
const poly &a,
const poly &b,
poly &c,
int var) {
832 bool both_mod_small=
false;
836 modq = (UWORD *)&c.modp;
838 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
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;
856 for (
int pow=maxpow; pow>=minpow; pow--) {
863 if (a[afirst+1+var] + b[blast+1+var] > pow) {
864 if (blast+b[blast] < b[0])
871 for (
int ai=afirst, bi=blast; ai<a[0] && bi>=1;) {
873 int thispow = AN.poly_num_vars==0 ? 0 : a[ai+1+var] + b[bi+1+var];
875 if (thispow == pow) {
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);
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);
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);
896 if (c.modp!=0) TakeNormalModulus((UWORD *)&c[ci+1+AN.poly_num_vars], &nc,
897 modq, nmodq, NOUNPACK);
901 bi -= ABS(b[bi-1]) + 2 + AN.poly_num_vars;
903 else if (thispow > pow)
906 bi -= ABS(b[bi-1]) + 2 + AN.poly_num_vars;
911 for (
int j=0; j<AN.poly_num_vars; j++)
913 if (AN.poly_num_vars > 0)
916 c[ci] = 2+AN.poly_num_vars+ABS(nc);
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]);
964 bool both_mod_small=
false;
968 modq = (UWORD *)&c.modp;
970 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
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;
986 for (
int i=0; i<AN.poly_num_vars; i++)
987 maxpowera[i] = maxpowerb[i] = 0;
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]);
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]);
997 for (
int i=0; i<AN.poly_num_vars; i++)
998 maxpower[i] = maxpowera[i] + maxpowerb[i];
1001 bool use_hash =
true;
1004 for (
int i=0; i<AN.poly_num_vars; i++) {
1005 if (nhash > POLY_MAX_HASH_SIZE / (maxpower[i]+1)) {
1010 nhash *= maxpower[i]+1;
1014 int nheap=a.number_of_terms();
1016 WantAddPointers(nheap+nhash);
1017 WORD **heap = AT.pWorkSpace + AT.pWorkPointer;
1019 for (
int ai=1, i=0; ai<a[0]; ai+=a[ai], i++) {
1020 heap[i] = (WORD *) NumberMalloc(
"poly::mul_heap");
1025 for (
int j=0; j<AN.poly_num_vars; j++)
1026 heap[i][4+j] = MAXPOSITIVE;
1029 WORD **hash = AT.pWorkSpace + AT.pWorkPointer + nheap;
1030 for (
int i=0; i<nhash; i++)
1040 pop_heap(BHEAD heap, nheap--);
1041 WORD *p = heap[nheap];
1045 if (use_hash) hash[p[2]] = NULL;
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]);
1058 ci = c.last_monomial_index();
1059 WORD nc = c[ci+c[ci]-1];
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)
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]);
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);
1079 if (c.modp!=0) TakeNormalModulus((UWORD *)&c[ci+1+AN.poly_num_vars], &nc,
1080 modq, nmodq, NOUNPACK);
1084 c[ci] = 2 + AN.poly_num_vars + ABS(nc);
1092 while (p[1] < b[0]) {
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];
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)
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]);
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]);
1116 if (c.modp!=0) TakeNormalModulus((UWORD *)&p[4+AN.poly_num_vars], &p[3],
1117 modq, nmodq, NOUNPACK);
1124 for (
int i=0; i<AN.poly_num_vars; i++)
1125 ID = (maxpower[i]+1)*ID + p[4+i];
1128 if (hash[ID] == NULL) {
1131 push_heap(BHEAD heap, ++nheap);
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)
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]);
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]);
1155 if (c.modp!=0) TakeNormalModulus((UWORD *)&h[4+AN.poly_num_vars], &h[3],
1156 modq, nmodq, NOUNPACK);
1163 push_heap(BHEAD heap, ++nheap);
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;
1196 if (a.is_zero() || b.is_zero()) { c[0]=1;
return; }
1198 c.check_memory(b[0]);
1199 c.termscopy(b.terms,0,b[0]);
1201 if (a.modp!=b.modp || a.modn!=b.modn) {
1203 c.setmod(a.modp,a.modn);
1208 c.check_memory(a[0]);
1209 c.termscopy(a.terms,0,a[0]);
1213 int na=a.number_of_terms();
1214 int nb=b.number_of_terms();
1216 if (na==1 || nb==1) {
1217 mul_one_term(a,b,c);
1224 if (vara!=-2 && varb!=-2 && vara==varb) {
1225 mul_univar(a,b,c,vara);
1241void poly::divmod_one_term (
const poly &a,
const poly &b,
poly &q,
poly &r,
bool only_divides) {
1243 POLY_GETIDENTITY(a);
1253 bool both_mod_small=
false;
1257 modq = (UWORD *)&q.modp;
1259 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
1260 both_mod_small=
true;
1266 ltbinv = NumberMalloc(
"poly::div_one_term");
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);
1274 GetLongModInverses(BHEAD (UWORD *)&b[2+AN.poly_num_vars], b[b[1]], modq, nmodq, ltbinv, &nltbinv, NULL, NULL);
1277 for (
int ai=1; ai<a[0]; ai+=a[ai]) {
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;
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);
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);
1309 MulLong((UWORD *)&a[ai+1+AN.poly_num_vars], a[ai+a[ai]-1],
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);
1323 r.termscopy(&a[ai+1+AN.poly_num_vars], ri+1+AN.poly_num_vars, ABS(nr));
1328 q[qi] = 2+AN.poly_num_vars+ABS(nq);
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]);
1344 if (only_divides) { r =
poly(BHEAD 1); ri=r[0];
break; }
1345 r[ri] = 2+AN.poly_num_vars+ABS(nr);
1354 if (q.modp!=0 || ltbinv != NULL) NumberFree(ltbinv,
"poly::div_one_term");
1374 POLY_GETIDENTITY(a);
1382 bool both_mod_small=
false;
1386 modq = (UWORD *)&q.modp;
1388 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
1389 both_mod_small=
true;
1394 ltbinv = NumberMalloc(
"poly::div_univar");
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);
1402 GetLongModInverses(BHEAD (UWORD *)&b[2+AN.poly_num_vars], b[b[1]], modq, nmodq, ltbinv, &nltbinv, NULL, NULL);
1407 UWORD *s = NumberMalloc(
"poly::div_univar");
1408 UWORD *t = NumberMalloc(
"poly::div_univar");
1411 int bpow = b[2+var];
1413 int ai=1, qi=1, ri=1;
1415 for (
int pow=a[2+var]; pow>=0; pow--) {
1421 while (ai<a[0] && a[ai+1+var] > pow)
1425 if (ai<a[0] && a[ai+1+var] == pow) {
1427 WCOPY(s, &a[ai+1+AN.poly_num_vars], ABS(ns));
1436 while (qj>1 && bi<b[0]) {
1438 qj -= 2 + AN.poly_num_vars + ABS(q[qj-1]);
1440 while (bi<b[0] && b[bi+1+var]+q[qj+1+var] > pow)
1443 if (bi<b[0] && b[bi+1+var]+q[qj+1+var] == pow) {
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);
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],
1457 AddLong(t,nt,s,ns,s,&ns);
1458 if (q.modp!=0) TakeNormalModulus((UWORD *)s,&ns,modq, nmodq, NOUNPACK);
1465 if (both_mod_small) {
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]);
1477 (UWORD *)&b[2+AN.poly_num_vars], b[b[1]],
1478 (UWORD *)&q[qi+1+AN.poly_num_vars], &ns, t, &nt);
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]);
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);
1505 for (
int i=0; i<AN.poly_num_vars; i++)
1507 q[qi+1+var] = pow-bpow;
1509 q[qi] = 2+AN.poly_num_vars+ABS(ns);
1515 if (only_divides) { r=
poly(BHEAD 1); ri=r[0];
break; }
1516 for (
int i=0; i<AN.poly_num_vars; i++)
1520 for (
int i=0; i<ABS(nt); i++)
1521 r[ri+1+AN.poly_num_vars+i] = t[i];
1523 r[ri] = 2+AN.poly_num_vars+ABS(nt);
1533 NumberFree(s,
"poly::div_univar");
1534 NumberFree(t,
"poly::div_univar");
1536 if (q.modp!=0 || ltbinv!=NULL) NumberFree(ltbinv,
"poly::div_univar");
1577 POLY_GETIDENTITY(a);
1587 LONG oldpWorkPointer = AT.pWorkPointer;
1589 bool both_mod_small=
false;
1593 modq = (UWORD *)&q.modp;
1595 if (a.modp>0 && b.modp>0 && a.modn==1 && b.modn==1)
1596 both_mod_small=
true;
1601 ltbinv = NumberMalloc(
"poly::div_heap-a");
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);
1609 GetLongModInverses(BHEAD (UWORD *)&b[2+AN.poly_num_vars], b[b[1]], modq, nmodq, ltbinv, &nltbinv, NULL, NULL);
1613 int nb=b.number_of_terms();
1614 WantAddPointers(nb);
1615 WORD **heap = AT.pWorkSpace + AT.pWorkPointer;
1616 AT.pWorkPointer += nb;
1618 for (
int i=0; i<nb; i++)
1619 heap[i] = (WORD *) NumberMalloc(
"poly::div_heap-b");
1624 WCOPY(&heap[0][3], &a[1], a[1]);
1625 heap[0][3] = a[a[1]];
1630 WORD *t = (WORD *) NumberMalloc(
"poly::div_heap-c");
1634 vector<pair<int,int> > insert;
1636 while (insert.size()>0 || nheap>0) {
1646 WORD *p = heap[nheap];
1649 if (insert.empty()) {
1651 this_insert =
false;
1653 pop_heap(BHEAD heap, nheap--);
1657 WCOPY(t, p, (5+ABS(p[3])+AN.poly_num_vars));
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)
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]);
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);
1686 p[0] = insert.back().first;
1687 p[1] = insert.back().second;
1696 if (p[0]==a[0])
break;
1697 WCOPY(&p[3], &a[p[0]], a[p[0]]);
1703 this_insert =
false;
1705 if (p[1]==qi) { s++;
break; }
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];
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)
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]);
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);
1739 swap (heap[nheap],p);
1740 push_heap(BHEAD heap, ++nheap);
1744 while (t[0]==-1 || (nheap>0 && monomial_compare(BHEAD heap[0]+3, t+3)==0));
1746 if (t[3] == 0)
continue;
1750 for (
int i=0; i<AN.poly_num_vars; i++)
1751 if (t[4+i] < b[2+i]) div=
false;
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]);
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)
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)
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]);
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);
1802 for (
int j=1; j<s; j++) {
1804 insert.push_back(make_pair(bi,qi));
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];
1816 r[ri] = 2+AN.poly_num_vars+ABS(nr);
1817 for (
int i=0; i<AN.poly_num_vars; i++)
1828 for (
int i=0; i<nb; i++)
1829 NumberFree(heap[i],
"poly::div_heap-b");
1831 NumberFree(t,
"poly::div_heap-c");
1833 if (q.modp!=0||ltbinv!=NULL) NumberFree(ltbinv,
"poly::div_heap-a");
1834 AT.pWorkPointer = oldpWorkPointer;
1854 q.modp = r.modp = a.modp;
1855 q.modn = r.modn = a.modn;
1863 MLOCK(ErrorMessageLock);
1864 MesPrint ((
char*)
"ERROR: division by zero polynomial");
1865 MUNLOCK(ErrorMessageLock);
1869 q.check_memory(a[0]);
1870 q.termscopy(a.terms,0,a[0]);
1875 if (b.is_monomial()) {
1876 divmod_one_term(a,b,q,r,only_divides);
1883 if (vara!=-2 && varb!=-2 && (vara==-1 || varb==-1 || vara==varb)) {
1898bool poly::divides (
const poly &a,
const poly &b) {
1899 POLY_GETIDENTITY(a);
1900 poly q(BHEAD 0), r(BHEAD 0);
1911void poly::div (
const poly &a,
const poly &b,
poly &c) {
1912 POLY_GETIDENTITY(a);
1923void poly::mod (
const poly &a,
const poly &b,
poly &c) {
1924 POLY_GETIDENTITY(a);
1935poly & poly::operator= (
const poly &a) {
1941 termscopy(a.terms,0,a[0]);
1953const poly poly::operator+ (
const poly &a)
const { POLY_GETIDENTITY(a);
poly b(BHEAD 0); add(*
this,a,b);
return b; }
1954const poly poly::operator- (
const poly &a)
const { POLY_GETIDENTITY(a);
poly b(BHEAD 0); sub(*
this,a,b);
return b; }
1955const poly poly::operator* (
const poly &a)
const { POLY_GETIDENTITY(a);
poly b(BHEAD 0);
mul(*
this,a,b);
return b; }
1956const poly poly::operator/ (
const poly &a)
const { POLY_GETIDENTITY(a);
poly b(BHEAD 0); div(*
this,a,b);
return b; }
1957const poly poly::operator% (
const poly &a)
const { POLY_GETIDENTITY(a);
poly b(BHEAD 0); mod(*
this,a,b);
return b; }
1960poly& poly::operator+= (
const poly &a) {
return *
this = *
this + a; }
1961poly& poly::operator-= (
const poly &a) {
return *
this = *
this - a; }
1962poly& poly::operator*= (
const poly &a) {
return *
this = *
this * a; }
1963poly& poly::operator/= (
const poly &a) {
return *
this = *
this / a; }
1964poly& poly::operator%= (
const poly &a) {
return *
this = *
this % a; }
1967bool poly::operator== (
const poly &a)
const {
1968 for (
int i=0; i<terms[0]; i++)
1969 if (terms[i] != a[i])
return 0;
1973bool poly::operator!= (
const poly &a)
const {
return !(*
this == a); }
1980int poly::number_of_terms ()
const {
1983 for (
int i=1; i<terms[0]; i+=terms[i])
1994int poly::first_variable ()
const {
1996 POLY_GETIDENTITY(*
this);
1998 int var = AN.poly_num_vars;
1999 for (
int j=0; j<var; j++)
2000 if (terms[2+j]>0)
return j;
2010const vector<int> poly::all_variables ()
const {
2012 POLY_GETIDENTITY(*
this);
2014 vector<bool> used(AN.poly_num_vars,
false);
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;
2021 for (
int i=0; i<AN.poly_num_vars; i++)
2022 if (used[i]) vars.push_back(i);
2033int poly::sign ()
const {
2034 if (terms[0]==1)
return 0;
2035 return terms[terms[1]] > 0 ? 1 : -1;
2044int poly::degree (
int x)
const {
2046 for (
int i=1; i<terms[0]; i+=terms[i])
2047 deg = MaX(deg, terms[i+1+x]);
2057int poly::total_degree ()
const {
2059 POLY_GETIDENTITY(*
this);
2062 for (
int i=1; i<terms[0]; i+=terms[i]) {
2064 for (
int j=0; j<AN.poly_num_vars; j++)
2065 deg += terms[i+1+j];
2066 tot_deg = MaX(deg, tot_deg);
2077const poly poly::integer_lcoeff ()
const {
2079 POLY_GETIDENTITY(*
this);
2085 res.termscopy(&terms[1],1,terms[1]);
2086 res[0] = res[1] + 1;
2087 for (
int i=0; i<AN.poly_num_vars; i++)
2099const poly poly::coefficient (
int x,
int n)
const {
2101 POLY_GETIDENTITY(*
this);
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;
2113 res[0] += res[res[0]];
2127const poly poly::lcoeff_univar (
int x)
const {
2128 return coefficient(x, degree(x));
2134const poly poly::lcoeff_multivar (
int x)
const {
2136 POLY_GETIDENTITY(*
this);
2138 poly res(BHEAD 0,modp,modn);
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;
2147 if (!same_powers)
break;
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;
2166const poly poly::derivative (
int x)
const {
2168 POLY_GETIDENTITY(*
this);
2173 for (
int i=1; i<terms[0]; i+=terms[i]) {
2175 int power = terms[i+1+x];
2179 b.termscopy(&terms[i], bi, terms[i]);
2182 WORD nb = b[bi+b[bi]-1];
2183 Product((UWORD *)&b[bi+1+AN.poly_num_vars], &nb, power);
2185 b[bi] = 2 + AN.poly_num_vars + ABS(nb);
2193 b.setmod(modp, modn);
2203bool poly::is_zero ()
const {
2204 return terms[0] == 1;
2213bool poly::is_one ()
const {
2215 POLY_GETIDENTITY(*
this);
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;
2233bool poly::is_integer ()
const {
2235 POLY_GETIDENTITY(*
this);
2237 if (terms[0] == 1)
return true;
2238 if (terms[0] != terms[1]+1)
return false;
2240 for (
int j=0; j<AN.poly_num_vars; j++)
2241 if (terms[2+j] != 0)
2253bool poly::is_monomial ()
const {
2254 return terms[0]>1 && terms[0]==terms[1]+1;
2280 POLY_GETIDENTITY(*
this);
2282 int num_terms=0, res=-1;
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;
2296 if (res == -1)
return -1;
2299 int deg = terms[2+res];
2300 if (2*num_terms < deg+1)
return -2;
2311const poly poly::simple_poly (PHEAD
int x,
int a,
int b,
int p,
int n) {
2313 poly tmp(BHEAD 0,p,n);
2316 tmp[idx++] = 3 + AN.poly_num_vars;
2317 for (
int i=0; i<AN.poly_num_vars; i++)
2318 tmp[idx++] = i==x ? 1 : 0;
2323 tmp[idx++] = 3 + AN.poly_num_vars;
2324 for (
int i=0; i<AN.poly_num_vars; i++) tmp[idx++] = 0;
2325 tmp[idx++] = ABS(a);
2326 tmp[idx++] = -SGN(a);
2331 if (b == 1)
return tmp;
2333 poly res(BHEAD 1,p,n);
2350const poly poly::simple_poly (PHEAD
int x,
const poly &a,
int b,
int p,
int n) {
2352 poly res(BHEAD 1,p,n);
2353 poly tmp(BHEAD 0,p,n);
2357 tmp[idx++] = 3 + AN.poly_num_vars;
2358 for (
int i=0; i<AN.poly_num_vars; i++)
2359 tmp[idx++] = i==x ? 1 : 0;
2364 tmp[idx++] = 2 + AN.poly_num_vars + ABS(a[a[0]-1]);
2365 for (
int i=0; i<AN.poly_num_vars; i++) tmp[idx++] = 0;
2366 for (
int i=0; i<ABS(a[a[0]-1]); i++)
2367 tmp[idx++] = a[2 + AN.poly_num_vars + i];
2368 tmp[idx++] = -a[a[0]-1];
2389void poly::get_variables (PHEAD vector<WORD *> es,
bool with_arghead,
bool sort_vars) {
2391 AN.poly_num_vars = 0;
2394 vector<int> degrees;
2395 map<int,int> var_to_idx;
2398 for (
int ei=0; ei<(int)es.size(); ei++) {
2402 if (*e == -SNUMBER) {
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);
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);
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]);
2427 degrees[var_to_idx[e[j]]] = MaX(degrees[var_to_idx[e[j]]], e[j+1]);
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]);
2448 if ( (AN.poly_num_vars+1)*
sizeof(WORD) > (
size_t)(AM.MaxTer) ) {
2450 AN.poly_vars = (WORD *)Malloc1((AN.poly_num_vars+1)*
sizeof(WORD),
"AN.poly_vars");
2451 AN.poly_vars_type = 1;
2454 AN.poly_vars = TermMalloc(
"AN.poly_vars");
2455 AN.poly_vars_type = 0;
2458 for (
int i=0; i<AN.poly_num_vars; i++)
2459 AN.poly_vars[i] = vars[i];
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]);
2473 sort(AN.poly_vars, AN.poly_vars+AN.poly_num_vars - var_to_idx.count(FACTORSYMBOL));
2483const poly poly::argument_to_poly (PHEAD WORD *e,
bool with_arghead,
bool sort_univar,
poly *denpoly) {
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;
2492 if (*e == -SNUMBER) {
2493 if (denpoly!=NULL) *denpoly =
poly(BHEAD 1);
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++)
2504 res[2+AN.poly_num_vars] = ABS(e[1]);
2505 res[3+AN.poly_num_vars] = SGN(e[1]);
2511 if (*e == -SYMBOL) {
2512 if (denpoly!=NULL) *denpoly =
poly(BHEAD 1);
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++)
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;
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");
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);
2541 if (denpoly!=NULL) *denpoly =
poly(BHEAD den, nden);
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];
2549 for (
int j=0; j<AN.poly_num_vars; j++)
2551 WCOPY(dum, &e[i+e[i]-ABS(nc)], ABS(nc));
2553 Mully(BHEAD dum, &nc, den, nden);
2554 res.termscopy((WORD *)dum, ri+1+AN.poly_num_vars, ABS(nc));
2555 res[ri] = ABS(nc) + AN.poly_num_vars + 2;
2556 res[ri+res[ri]-1] = nc;
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];
2567 if (sort_univar || AN.poly_num_vars>1)
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");
2584void poly::poly_to_argument (
const poly &a, WORD *res,
bool with_arghead) {
2586 POLY_GETIDENTITY(a);
2601 res[1] = AN.poly_num_vars>1 ? DIRTYFLAG : 0;
2602 for (
int i=2; i<ARGHEAD; i++)
2606 int L = with_arghead ? ARGHEAD : 0;
2608 for (
int i=1; i!=a[0]; i+=a[i]) {
2614 for (
int j=0; j<AN.poly_num_vars; j++)
2621 res[L+1+res[L+2]++] = AN.poly_vars[j];
2622 res[L+1+res[L+2]++] = a[i+1+j];
2625 if (!first) res[L] += res[L+2];
2627 WORD nc = a[i+a[i]-1];
2628 WCOPY(&res[L+res[L]], &a[i+a[i]-1-ABS(nc)], ABS(nc));
2631 memset(&res[L+res[L]], 0, ABS(nc)*
sizeof(WORD));
2634 res[L+res[L]] = SGN(nc) * (2*ABS(nc)+1);
2656void poly::poly_to_argument_with_den (
const poly &a, WORD nden,
const UWORD *den, WORD *res,
bool with_arghead) {
2658 POLY_GETIDENTITY(a);
2673 res[1] = AN.poly_num_vars>1 ? DIRTYFLAG : 0;
2674 for (
int i=2; i<ARGHEAD; i++)
2679 UWORD *den1 = (UWORD *)NumberMalloc(
"poly_to_argument_with_den");
2681 int L = with_arghead ? ARGHEAD : 0;
2683 for (
int i=1; i!=a[0]; i+=a[i]) {
2689 for (
int j=0; j<AN.poly_num_vars; j++)
2696 res[L+1+res[L+2]++] = AN.poly_vars[j];
2697 res[L+1+res[L+2]++] = a[i+1+j];
2700 if (!first) res[L] += res[L+2];
2703 WORD nc = a[i+a[i]-1];
2704 WCOPY(&res[L+res[L]], &a[i+a[i]-1-ABS(nc)], ABS(nc));
2708 WCOPY(den1, den, ABS(nden));
2710 if (nden != 1 || den[0] != 1) {
2712 Simplify(BHEAD (UWORD *)&res[L+res[L]], &nc, den1, &nden1);
2715 Pack((UWORD *)&res[L+res[L]], &nc, den1, nden1);
2716 res[L] += 2*ABS(nc)+1;
2717 res[L+res[L]-1] = SGN(nc)*(2*ABS(nc)+1);
2721 NumberFree(den1,
"poly_to_argument_with_den");
2739int poly::size_of_form_notation ()
const {
2741 POLY_GETIDENTITY(*
this);
2744 if (terms[0]==1)
return 0;
2748 for (
int i=1; i!=terms[0]; i+=terms[i]) {
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;
2768int poly::size_of_form_notation_with_den (WORD nden)
const {
2770 POLY_GETIDENTITY(*
this);
2773 if (terms[0]==1)
return 0;
2778 for (
int i=1; i!=terms[0]; i+=terms[i]) {
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;
2796const vector<WORD> poly::to_coefficient_list (
const poly &a) {
2798 POLY_GETIDENTITY(a);
2800 if (a.is_zero())
return vector<WORD>();
2802 int x = a.first_variable();
2803 if (x == AN.poly_num_vars) x=0;
2805 vector<WORD> res(1+a[2+x],0);
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;
2819const vector<WORD> poly::coefficient_list_divmod (
const vector<WORD> &a,
const vector<WORD> &b, WORD p,
int divmod) {
2821 int bsize = (int)b.size();
2822 while (b[bsize-1]==0) bsize--;
2827 vector<WORD> q(a.size(),0);
2830 while ((
int)r.size() >= bsize) {
2831 LONG
mul = ((LONG)inv * r.back()) % p;
2832 int off = r.size()-bsize;
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)
2841 while (q.size()>0 && q.back()==0)
2846 while (r.size()>0 && r.back()==0)
2858const poly poly::from_coefficient_list (PHEAD
const vector<WORD> &a,
int x, WORD p) {
2863 for (
int i=(
int)a.size()-1; i>=0; i--)
2865 res.check_memory(ri);
2866 res[ri] = AN.poly_num_vars+3;
2867 for (
int j=0; j<AN.poly_num_vars; j++)
2870 res[ri+1+AN.poly_num_vars] = ABS(a[i]);
2871 res[ri+1+AN.poly_num_vars+1] = SGN(a[i]);
static void divmod(const poly &, const poly &, poly &, poly &, bool only_divides)
static void mul_heap(const poly &, const poly &, poly &)
int is_dense_univariate() const
static void divmod_univar(const poly &, const poly &, poly &, poly &, int, bool)
static void mul(const poly &, const poly &, poly &)
static void divmod_heap(const poly &, const poly &, poly &, poly &, bool, bool, bool &)
VOID RaisPowCached(PHEAD WORD, WORD, UWORD **, WORD *)
int GetModInverses(WORD, WORD, WORD *, WORD *)