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