FORM  4.3
polywrap.cc
Go to the documentation of this file.
1 
8 /* #[ License : */
9 /*
10  * Copyright (C) 1984-2022 J.A.M. Vermaseren
11  * When using this file you are requested to refer to the publication
12  * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
13  * This is considered a matter of courtesy as the development was paid
14  * for by FOM the Dutch physics granting agency and we would like to
15  * be able to track its scientific use to convince FOM of its value
16  * for the community.
17  *
18  * This file is part of FORM.
19  *
20  * FORM is free software: you can redistribute it and/or modify it under the
21  * terms of the GNU General Public License as published by the Free Software
22  * Foundation, either version 3 of the License, or (at your option) any later
23  * version.
24  *
25  * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
26  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
27  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
28  * details.
29  *
30  * You should have received a copy of the GNU General Public License along
31  * with FORM. If not, see <http://www.gnu.org/licenses/>.
32  */
33 /* #] License : */
34 
35 #include "poly.h"
36 #include "polygcd.h"
37 #include "polyfact.h"
38 
39 #include <iostream>
40 #include <vector>
41 #include <map>
42 #include <climits>
43 #include <cassert>
44 
45 //#define DEBUG
46 
47 #ifdef DEBUG
48 #include "mytime.h"
49 #endif
50 
51 // denompower is increased by this factor when divmod_heap fails
52 const int POLYWRAP_DENOMPOWER_INCREASE_FACTOR = 2;
53 
54 using namespace std;
55 
56 /*
57  #[ poly_determine_modulus :
58 */
59 
79 WORD poly_determine_modulus (PHEAD bool multi_error, bool is_fun_arg, string message) {
80 
81  if (AC.ncmod==0) return 0;
82 
83  if (!is_fun_arg || (AC.modmode & ALSOFUNARGS)) {
84 
85  if (ABS(AC.ncmod)>1) {
86  MLOCK(ErrorMessageLock);
87  MesPrint ((char*)"ERROR: %s with modulus > WORDSIZE not implemented",message.c_str());
88  MUNLOCK(ErrorMessageLock);
89  Terminate(-1);
90  }
91 
92  if (multi_error && AN.poly_num_vars>1) {
93  MLOCK(ErrorMessageLock);
94  MesPrint ((char*)"ERROR: multivariate %s with modulus not implemented",message.c_str());
95  MUNLOCK(ErrorMessageLock);
96  Terminate(-1);
97  }
98 
99  return *AC.cmod;
100  }
101 
102  AN.ncmod = 0;
103  return 0;
104 }
105 
106 /*
107  #] poly_determine_modulus :
108  #[ poly_gcd :
109 */
110 
124 WORD *poly_gcd(PHEAD WORD *a, WORD *b, WORD fit) {
125 
126 #ifdef DEBUG
127  cout << "*** [" << thetime() << "] CALL : poly_gcd" << endl;
128 #endif
129 
130 //
131 //MesPrint("Calling poly_gcd with:");
132 //{
133 // WORD *at = a;
134 // MesPrint(" a:");
135 // while ( *at ) {
136 // MesPrint(" %a",*at,at);
137 // at += *at;
138 // }
139 // MesPrint(" b:");
140 // at = b;
141 // while ( *at ) {
142 // MesPrint(" %a",*at,at);
143 // at += *at;
144 // }
145 //}
146 
147  // Extract variables
148  vector<WORD *> e;
149  e.push_back(a);
150  e.push_back(b);
151  poly::get_variables(BHEAD e, false, true);
152 
153  // Check for modulus calculus
154  WORD modp=poly_determine_modulus(BHEAD true, true, "polynomial GCD");
155 
156  // Convert to polynomials
157  poly pa(poly::argument_to_poly(BHEAD a, false, true), modp, 1);
158  poly pb(poly::argument_to_poly(BHEAD b, false, true), modp, 1);
159 
160  // Calculate gcd
161  poly gcd(polygcd::gcd(pa,pb));
162 
163  // Allocate new memory and convert to Form notation
164  int newsize = (gcd.size_of_form_notation()+1);
165  WORD *res;
166  if ( fit ) {
167  if ( newsize*sizeof(WORD) >= (size_t)(AM.MaxTer) ) {
168  MLOCK(ErrorMessageLock);
169  MesPrint("poly_gcd: Term too complex. Maybe increasing MaxTermSize can help");
170  MUNLOCK(ErrorMessageLock);
171  Terminate(-1);
172  }
173  res = TermMalloc("poly_gcd");
174  }
175  else {
176  res = (WORD *)Malloc1(newsize*sizeof(WORD), "poly_gcd");
177  }
178  poly::poly_to_argument(gcd, res, false);
179 
180  poly_free_poly_vars(BHEAD "AN.poly_vars_qcd");
181 
182  // reset modulo calculation
183  AN.ncmod = AC.ncmod;
184  return res;
185 }
186 
187 /*
188  #] poly_gcd :
189  #[ poly_divmod :
190 
191  if fit == 1 the answer must fit inside a term.
192 */
193 
194 WORD *poly_divmod(PHEAD WORD *a, WORD *b, int divmod, WORD fit) {
195 
196 #ifdef DEBUG
197  cout << "*** [" << thetime() << "] CALL : poly_divmod" << endl;
198 #endif
199 
200  // check for modulus calculus
201  WORD modp=poly_determine_modulus(BHEAD false, true, "polynomial division");
202 
203  // get variables
204  vector<WORD *> e;
205  e.push_back(a);
206  e.push_back(b);
207  poly::get_variables(BHEAD e, false, false);
208 
209  // add extra variables to keep track of denominators
210  const int DENOMSYMBOL = MAXPOSITIVE;
211 
212 // WORD *new_poly_vars = (WORD *)Malloc1((AN.poly_num_vars+1)*sizeof(WORD), "AN.poly_vars");
213 // WCOPY(new_poly_vars, AN.poly_vars, AN.poly_num_vars);
214 // new_poly_vars[AN.poly_num_vars] = DENOMSYMBOL;
215 // if (AN.poly_num_vars > 0)
216 // M_free(AN.poly_vars, "AN.poly_vars");
217 // AN.poly_num_vars++;
218 // AN.poly_vars = new_poly_vars;
219  AN.poly_vars[AN.poly_num_vars++] = DENOMSYMBOL;
220 
221  // convert to polynomials
222  poly dena(BHEAD 0);
223  poly denb(BHEAD 0);
224  poly pa(poly::argument_to_poly(BHEAD a, false, true, &dena), modp, 1);
225  poly pb(poly::argument_to_poly(BHEAD b, false, true, &denb), modp, 1);
226 
227  // remove contents
228  poly numres(polygcd::integer_content(pa));
229  poly denres(polygcd::integer_content(pb));
230  pa /= numres;
231  pb /= denres;
232 
233  if (divmod==0) {
234  numres *= denb;
235  denres *= dena;
236  }
237  else {
238  denres = dena;
239  }
240 
241  poly gcdres(polygcd::integer_gcd(numres,denres));
242  numres /= gcdres;
243  denres /= gcdres;
244 
245  // determine lcoeff(b)
246  poly lcoeffb(pb.integer_lcoeff());
247 
248  poly pres(BHEAD 0);
249 
250  if (!lcoeffb.is_one()) {
251 
252  if (AN.poly_num_vars > 2) {
253  // the original polynomial is multivariate (one dummy variable has
254  // been added), so it is not trivial to determine which power of
255  // lcoeff(b) can be in the answer
256 
257  int denompower = 0, prevdenompower = 0;
258  poly pq(BHEAD 0);
259  poly pr(BHEAD 0);
260 
261  // try denompower = 0, if this fails increase denompower until division succeeds
262  bool div_fail = true;
263  do
264  {
265  if(denompower < prevdenompower)
266  {
267  // denompower increased beyond INT_MAX
268  MLOCK(ErrorMessageLock);
269  MesPrint ((char*)"ERROR: pseudo-division failed in poly_divmod (denompower > INT_MAX)");
270  MUNLOCK(ErrorMessageLock);
271  Terminate(1);
272  }
273 
274  if(denompower != 0)
275  {
276  // multiply a by lcoeffb^(denompower-prevdenompower)
277  WORD n = lcoeffb[lcoeffb[1]];
278  RaisPow(BHEAD (UWORD *)&lcoeffb[2+AN.poly_num_vars], &n, denompower-prevdenompower);
279  lcoeffb[1] = 2 + AN.poly_num_vars + ABS(n);
280  lcoeffb[0] = 1 + lcoeffb[1];
281  lcoeffb[lcoeffb[1]] = n;
282 
283  pa *= lcoeffb;
284  denres *= lcoeffb;
285  }
286 
287  // try division
288  poly ppow(BHEAD 0);
289  poly::divmod_heap(pa,pb,pq,pr,false,true,div_fail); // sets div_fail
290 
291  // increase denompower for next iteration
292  prevdenompower = denompower;
293  denompower = (denompower==0 ? 1 : denompower*POLYWRAP_DENOMPOWER_INCREASE_FACTOR+1 ); // generates 2^n-1 when POLYWRAP_DENOMPOWER_INCREASE_FACTOR = 2
294  }
295  while(div_fail);
296 
297  pres = (divmod==0 ? pq : pr);
298 
299  }
300  else {
301  // one variable, so the power is the difference of the degrees
302 
303  int denompower = MaX(0, pa.degree(0) - pb.degree(0) + 1);
304 
305  // multiply a by that power
306  WORD n = lcoeffb[lcoeffb[1]];
307  RaisPow(BHEAD (UWORD *)&lcoeffb[2+AN.poly_num_vars], &n, denompower);
308  lcoeffb[1] = 2 + AN.poly_num_vars + ABS(n);
309  lcoeffb[0] = 1 + lcoeffb[1];
310  lcoeffb[lcoeffb[1]] = n;
311 
312  pa *= lcoeffb;
313  denres *= lcoeffb;
314 
315  pres = (divmod==0 ? pa/pb : pa%pb);
316 
317  }
318 
319  }
320  else {
321 
322  pres = (divmod==0 ? pa/pb : pa%pb);
323 
324  }
325 
326  // convert to Form notation
327  // NOTE: this part can be rewritten with poly::size_of_form_notation_with_den()
328  // and poly::poly_to_argument_with_den().
329  WORD *res;
330 
331  // special case: a=0
332  if (pres[0]==1) {
333  if ( fit ) {
334  res = TermMalloc("poly_divmod");
335  }
336  else {
337  res = (WORD *)Malloc1(sizeof(WORD), "poly_divmod");
338  }
339  res[0] = 0;
340  }
341  else {
342  pres *= numres;
343 
344  WORD nden;
345  UWORD *den = (UWORD *)NumberMalloc("poly_divmod");
346 
347  // allocate the memory; note that this overestimates the size,
348  // since the estimated denominators are too large
349  int ressize = pres.size_of_form_notation() + pres.number_of_terms()*2*ABS(denres[denres[1]]) + 1;
350 
351  if ( fit ) {
352  if ( ressize*sizeof(WORD) > (size_t)(AM.MaxTer) ) {
353  MLOCK(ErrorMessageLock);
354  MesPrint("poly_divmod: Term too complex. Maybe increasing MaxTermSize can help");
355  MUNLOCK(ErrorMessageLock);
356  Terminate(-1);
357  }
358  res = TermMalloc("poly_divmod");
359  }
360  else {
361  res = (WORD *)Malloc1(ressize*sizeof(WORD), "poly_divmod");
362  }
363  int L=0;
364 
365  for (int i=1; i!=pres[0]; i+=pres[i]) {
366 
367  res[L]=1; // length
368  bool first = true;
369 
370  for (int j=0; j<AN.poly_num_vars; j++)
371  if (pres[i+1+j] > 0) {
372  if (first) {
373  first = false;
374  res[L+1] = 1; // symbols
375  res[L+2] = 2; // length
376  }
377  res[L+1+res[L+2]++] = AN.poly_vars[j]; // symbol
378  res[L+1+res[L+2]++] = pres[i+1+j]; // power
379  }
380 
381  if (!first) res[L] += res[L+2]; // fix length
382 
383  // numerator
384  WORD nnum = pres[i+pres[i]-1];
385  WCOPY(&res[L+res[L]], &pres[i+pres[i]-1-ABS(nnum)], ABS(nnum));
386 
387  // calculate denominator
388  nden = denres[denres[1]];
389  WCOPY(den, &denres[2+AN.poly_num_vars], ABS(nden));
390 
391  if (nden!=1 || den[0]!=1)
392  Simplify(BHEAD (UWORD *)&res[L+res[L]], &nnum, den, &nden); // gcd(num,den)
393  Pack((UWORD *)&res[L+res[L]], &nnum, den, nden); // format
394  res[L] += 2*ABS(nnum)+1; // fix length
395  res[L+res[L]-1] = SGN(nnum)*(2*ABS(nnum)+1); // length of coefficient
396  L += res[L]; // fix length
397  }
398 
399  res[L] = 0;
400 
401  NumberFree(den,"poly_divmod");
402  }
403 
404  // clean up
405  poly_free_poly_vars(BHEAD "AN.poly_vars_divmod");
406 
407  // reset modulo calculation
408  AN.ncmod = AC.ncmod;
409 
410  return res;
411 }
412 
413 /*
414  #] poly_divmod :
415  #[ poly_div :
416 
417  Routine divides the expression in arg1 by the expression in arg2.
418  We did not take out special cases.
419  The arguments are zero terminated sequences of term(s).
420  The action is to divide arg1 by arg2: [arg1/arg2].
421  The answer should be a buffer (allocated by Malloc1) with a zero
422  terminated sequence of terms (or just zero).
423 */
424 WORD *poly_div(PHEAD WORD *a, WORD *b, WORD fit) {
425 
426 #ifdef DEBUG
427  cout << "*** [" << thetime() << "] CALL : poly_div" << endl;
428 #endif
429 
430  return poly_divmod(BHEAD a, b, 0, fit);
431 }
432 
433 /*
434  #] poly_div :
435  #[ poly_rem :
436 
437  Routine divides the expression in arg1 by the expression in arg2
438  and takes the remainder.
439  We did not take out special cases.
440  The arguments are zero terminated sequences of term(s).
441  The action is to divide arg1 by arg2 and take the remainder: [arg1%arg2].
442  The answer should be a buffer (allocated by Malloc1) with a zero
443  terminated sequence of terms (or just zero).
444 */
445 WORD *poly_rem(PHEAD WORD *a, WORD *b, WORD fit) {
446 
447 #ifdef DEBUG
448  cout << "*** [" << thetime() << "] CALL : poly_rem" << endl;
449 #endif
450 
451  return poly_divmod(BHEAD a, b, 1, fit);
452 }
453 
454 /*
455  #] poly_rem :
456  #[ poly_ratfun_read :
457 */
458 
471 void poly_ratfun_read (WORD *a, poly &num, poly &den) {
472 
473 #ifdef DEBUG
474  cout << "*** [" << thetime() << "] CALL : poly_ratfun_read" << endl;
475 #endif
476 
477  POLY_GETIDENTITY(num);
478 
479  int modp = num.modp;
480 
481  WORD *astop = a+a[1];
482 
483  bool clean = (a[2] & MUSTCLEANPRF) == 0;
484 
485  a += FUNHEAD;
486  if (a >= astop) {
487  MLOCK(ErrorMessageLock);
488  MesPrint ((char*)"ERROR: PolyRatFun cannot have zero arguments");
489  MUNLOCK(ErrorMessageLock);
490  Terminate(-1);
491  }
492 
493  poly den_num(BHEAD 1),den_den(BHEAD 1);
494 
495  num = poly::argument_to_poly(BHEAD a, true, !clean, &den_num);
496  num.setmod(modp,1);
497  NEXTARG(a);
498 
499  if (a < astop) {
500  den = poly::argument_to_poly(BHEAD a, true, !clean, &den_den);
501  den.setmod(modp,1);
502  NEXTARG(a);
503  }
504  else {
505  den = poly(BHEAD 1, modp, 1);
506  }
507 
508  if (a < astop) {
509  MLOCK(ErrorMessageLock);
510  MesPrint ((char*)"ERROR: PolyRatFun cannot have more than two arguments");
511  MUNLOCK(ErrorMessageLock);
512  Terminate(-1);
513  }
514 
515  if (!clean) {
516  vector<WORD> minpower(AN.poly_num_vars, MAXPOSITIVE);
517 
518  for (int i=1; i<num[0]; i+=num[i])
519  for (int j=0; j<AN.poly_num_vars; j++)
520  minpower[j] = MiN(minpower[j], num[i+1+j]);
521  for (int i=1; i<den[0]; i+=den[i])
522  for (int j=0; j<AN.poly_num_vars; j++)
523  minpower[j] = MiN(minpower[j], den[i+1+j]);
524 
525  for (int i=1; i<num[0]; i+=num[i])
526  for (int j=0; j<AN.poly_num_vars; j++)
527  num[i+1+j] -= minpower[j];
528  for (int i=1; i<den[0]; i+=den[i])
529  for (int j=0; j<AN.poly_num_vars; j++)
530  den[i+1+j] -= minpower[j];
531 
532  num *= den_den;
533  den *= den_num;
534 
535  poly gcd = polygcd::gcd(num,den);
536  num /= gcd;
537  den /= gcd;
538  }
539 }
540 
541 /*
542  #] poly_ratfun_read :
543  #[ poly_sort :
544 */
545 
557 void poly_sort(PHEAD WORD *a) {
558 
559 #ifdef DEBUG
560  cout << "*** [" << thetime() << "] CALL : poly_sort" << endl;
561 #endif
562  if (NewSort(BHEAD0)) { Terminate(-1); }
563  AR.CompareRoutine = &CompareSymbols;
564 
565  for (int i=ARGHEAD; i<a[0]; i+=a[i]) {
566  if (SymbolNormalize(a+i)<0 || StoreTerm(BHEAD a+i)) {
567  AR.CompareRoutine = &Compare1;
568  LowerSortLevel();
569  Terminate(-1);
570  }
571  }
572 
573  if (EndSort(BHEAD a+ARGHEAD,1) < 0) {
574  AR.CompareRoutine = &Compare1;
575  Terminate(-1);
576  }
577 
578  AR.CompareRoutine = &Compare1;
579  a[1] = 0; // set dirty flag to zero
580 }
581 
582 /*
583  #] poly_sort :
584  #[ poly_ratfun_add :
585 */
586 
600 WORD *poly_ratfun_add (PHEAD WORD *t1, WORD *t2) {
601 
602  if ( AR.PolyFunExp == 1 ) return PolyRatFunSpecial(BHEAD t1, t2);
603 
604 #ifdef DEBUG
605  cout << "*** [" << thetime() << "] CALL : poly_ratfun_add" << endl;
606 #endif
607 
608  WORD *oldworkpointer = AT.WorkPointer;
609 
610  // Extract variables
611  vector<WORD *> e;
612 
613  for (WORD *t=t1+FUNHEAD; t<t1+t1[1];) {
614  e.push_back(t);
615  NEXTARG(t);
616  }
617  for (WORD *t=t2+FUNHEAD; t<t2+t2[1];) {
618  e.push_back(t);
619  NEXTARG(t);
620  }
621 
622  poly::get_variables(BHEAD e, true, true);
623 
624  // Check for modulus calculus
625  WORD modp=poly_determine_modulus(BHEAD true, true, "PolyRatFun");
626 
627  // Find numerators / denominators
628  poly num1(BHEAD 0,modp,1), den1(BHEAD 0,modp,1), num2(BHEAD 0,modp,1), den2(BHEAD 0,modp,1);
629 
630  poly_ratfun_read(t1, num1, den1);
631  poly_ratfun_read(t2, num2, den2);
632 
633  poly num(BHEAD 0),den(BHEAD 0),gcd(BHEAD 0);
634 
635  // Calculate result
636  if (den1 != den2) {
637  gcd = polygcd::gcd(den1,den2);
638 #ifdef OLDADDITION
639  num = num1*(den2/gcd) + num2*(den1/gcd);
640  den = (den1/gcd)*den2;
641  gcd = polygcd::gcd(num,den);
642 #else
643  den = den1/gcd;
644  num = num1*(den2/gcd) + num2*den;
645  den = den*den2;
646  gcd = polygcd::gcd(num,gcd);
647 #endif
648  }
649  else {
650  num = num1 + num2;
651  den = den1;
652  gcd = polygcd::gcd(num,den);
653  }
654 
655  num /= gcd;
656  den /= gcd;
657 
658  // Fix sign
659  if (den.sign() == -1) { num*=poly(BHEAD -1); den*=poly(BHEAD -1); }
660 
661  // Check size
662  if (num.size_of_form_notation() + den.size_of_form_notation() + 3 >= AM.MaxTer/(int)sizeof(WORD)) {
663  MLOCK(ErrorMessageLock);
664  MesPrint ("ERROR: PolyRatFun doesn't fit in a term");
665  MesPrint ("(1) num size = %d, den size = %d, MaxTer = %d",num.size_of_form_notation(),
666  den.size_of_form_notation(),AM.MaxTer);
667  MUNLOCK(ErrorMessageLock);
668  Terminate(-1);
669  }
670 
671  // Format result in Form notation
672  WORD *t = oldworkpointer;
673 
674  *t++ = AR.PolyFun; // function
675  *t++ = 0; // length (to be determined)
676 // *t++ &= ~MUSTCLEANPRF; // clean polyratfun
677  *t++ = 0;
678  FILLFUN3(t); // header
679  poly::poly_to_argument(num,t, true); // argument 1 (numerator)
680  if (*t>0 && t[1]==DIRTYFLAG) // to Form order
681  poly_sort(BHEAD t);
682  t += (*t>0 ? *t : 2);
683  poly::poly_to_argument(den,t, true); // argument 2 (denominator)
684  if (*t>0 && t[1]==DIRTYFLAG) // to Form order
685  poly_sort(BHEAD t);
686  t += (*t>0 ? *t : 2);
687 
688  oldworkpointer[1] = t - oldworkpointer; // length
689  AT.WorkPointer = t;
690 
691  poly_free_poly_vars(BHEAD "AN.poly_vars_ratfun_add");
692 
693  // reset modulo calculation
694  AN.ncmod = AC.ncmod;
695 
696  return oldworkpointer;
697 }
698 
699 /*
700  #] poly_ratfun_add :
701  #[ poly_ratfun_normalize :
702 */
703 
719 int poly_ratfun_normalize (PHEAD WORD *term) {
720 
721 #ifdef DEBUG
722  cout << "*** [" << thetime() << "] CALL : poly_ratfun_normalize" << endl;
723 #endif
724 
725  // Strip coefficient
726  WORD *tstop = term + *term;
727  int ncoeff = tstop[-1];
728  tstop -= ABS(ncoeff);
729 
730  // if only one clean polyratfun, return immediately
731  int num_polyratfun = 0;
732 
733  for (WORD *t=term+1; t<tstop; t+=t[1])
734  if (*t == AR.PolyFun) {
735  num_polyratfun++;
736  if ((t[2] & MUSTCLEANPRF) != 0)
737  num_polyratfun = INT_MAX;
738  if (num_polyratfun > 1) break;
739  }
740 
741  if (num_polyratfun <= 1) return 0;
742 
743  WORD oldsorttype = AR.SortType;
744  AR.SortType = SORTHIGHFIRST;
745 
746 /*
747  When there are polyratfun's with only one variable: rename them
748  temporarily to TMPPOLYFUN.
749 */
750  for (WORD *t=term+1; t<tstop; t+=t[1]) {
751  if (*t == AR.PolyFun && (t[1] == FUNHEAD+t[FUNHEAD]
752  || t[1] == FUNHEAD+2 ) ) { *t = TMPPOLYFUN; }
753  }
754 
755 
756  // Extract all variables in the polyfuns
757  vector<WORD *> e;
758 
759  for (WORD *t=term+1; t<tstop; t+=t[1]) {
760  if (*t == AR.PolyFun)
761  for (WORD *t2 = t+FUNHEAD; t2<t+t[1];) {
762  e.push_back(t2);
763  NEXTARG(t2);
764  }
765  }
766  poly::get_variables(BHEAD e, true, true);
767 
768  // Check for modulus calculus
769  WORD modp=poly_determine_modulus(BHEAD true, true, "PolyRatFun");
770 
771  // Accumulate total denominator/numerator and copy the remaining terms
772  // We start with 'trivial' polynomials
773  poly num1(BHEAD (UWORD *)tstop, ncoeff/2, modp, 1);
774  poly den1(BHEAD (UWORD *)tstop+ABS(ncoeff/2), ABS(ncoeff)/2, modp, 1);
775 
776  WORD *s = term+1;
777 
778  for (WORD *t=term+1; t<tstop;)
779  if (*t == AR.PolyFun) {
780 
781  poly num2(BHEAD 0,modp,1);
782  poly den2(BHEAD 0,modp,1);
783  poly_ratfun_read(t,num2,den2);
784  if ((t[2] & MUSTCLEANPRF) != 0) { // first normalize
785  poly gcd1(polygcd::gcd(num2,den2));
786  num2 = num2/gcd1;
787  den2 = den2/gcd1;
788  }
789  t += t[1];
790  poly gcd1(polygcd::gcd(num1,den2));
791  poly gcd2(polygcd::gcd(num2,den1));
792 
793  num1 = (num1 / gcd1) * (num2 / gcd2);
794  den1 = (den1 / gcd2) * (den2 / gcd1);
795  }
796  else {
797  int i = t[1];
798  if ( s != t ) { NCOPY(s,t,i) }
799  else { t += i; s += i; }
800  }
801 
802  // Fix sign
803  if (den1.sign() == -1) { num1*=poly(BHEAD -1); den1*=poly(BHEAD -1); }
804 
805  // Check size
806  if (num1.size_of_form_notation() + den1.size_of_form_notation() + 3 >= AM.MaxTer/(int)sizeof(WORD)) {
807  MLOCK(ErrorMessageLock);
808  MesPrint ("ERROR: PolyRatFun doesn't fit in a term");
809  MesPrint ("(2) num size = %d, den size = %d, MaxTer = %d",num1.size_of_form_notation(),
810  den1.size_of_form_notation(),AM.MaxTer);
811  MUNLOCK(ErrorMessageLock);
812  Terminate(-1);
813  }
814 
815  // Format result in Form notation
816  WORD *t = s;
817  *t++ = AR.PolyFun; // function
818  *t++ = 0; // size (to be determined)
819  *t++ &= ~MUSTCLEANPRF; // clean polyratfun
820  FILLFUN3(t); // header
821  poly::poly_to_argument(num1,t,true); // argument 1 (numerator)
822  if (*t>0 && t[1]==DIRTYFLAG) // to Form order
823  poly_sort(BHEAD t);
824  t += (*t>0 ? *t : 2);
825  poly::poly_to_argument(den1,t,true); // argument 2 (denominator)
826  if (*t>0 && t[1]==DIRTYFLAG) // to Form order
827  poly_sort(BHEAD t);
828  t += (*t>0 ? *t : 2);
829 
830  s[1] = t - s; // function length
831 
832  *t++ = 1; // term coefficient
833  *t++ = 1;
834  *t++ = 3;
835 
836  term[0] = t-term; // term length
837 
838  poly_free_poly_vars(BHEAD "AN.poly_vars_ratfun_normalize");
839 
840  // reset modulo calculation
841  AN.ncmod = AC.ncmod;
842 
843  tstop = term + *term; tstop -= ABS(tstop[-1]);
844  for (WORD *t=term+1; t<tstop; t+=t[1]) {
845  if (*t == TMPPOLYFUN ) *t = AR.PolyFun;
846  }
847 
848  AR.SortType = oldsorttype;
849  return 0;
850 }
851 
852 /*
853  #] poly_ratfun_normalize :
854  #[ poly_fix_minus_signs :
855 */
856 
857 void poly_fix_minus_signs (factorized_poly &a) {
858 
859  if ( a.factor.empty() ) return;
860 
861  POLY_GETIDENTITY(a.factor[0]);
862 
863  int overall_sign = +1;
864 
865  // find term with maximum power of highest symbol
866  for (int i=0; i<(int)a.factor.size(); i++) {
867 
868  int maxvar = -1;
869  int maxpow = -1;
870  int sign = +1;
871 
872  WORD *tstop = a.factor[i].terms; tstop += *tstop;
873  for (WORD *t=a.factor[i].terms+1; t<tstop; t+=*t)
874  for (int j=0; j<AN.poly_num_vars; j++) {
875  int var = AN.poly_vars[j];
876  int pow = t[1+j];
877  if (pow>0 && (var>maxvar || (var==maxvar && pow>maxpow))) {
878  maxvar = var;
879  maxpow = pow;
880  sign = SGN(*(t+*t-1));
881  }
882  }
883 
884  // if negative coefficient, multiply by -1
885  if (sign==-1) {
886  a.factor[i] *= poly(BHEAD sign);
887  if (a.power[i] % 2 == 1) overall_sign*=-1;
888  }
889  }
890 
891  // if overall minus sign
892  if (overall_sign == -1) {
893  // look at constant factor and multiply by -1
894  for (int i=0; i<(int)a.factor.size(); i++)
895  if (a.factor[i].is_integer()) {
896  a.factor[i] *= poly(BHEAD -1);
897  return;
898  }
899 
900  // otherwise, add a factor of -1
901  a.add_factor(poly(BHEAD -1), 1);
902  }
903 }
904 
905 /*
906  #] poly_fix_minus_signs :
907  #[ poly_factorize :
908 */
909 
922 WORD *poly_factorize (PHEAD WORD *argin, WORD *argout, bool with_arghead, bool is_fun_arg) {
923 
924 #ifdef DEBUG
925  cout << "*** [" << thetime() << "] CALL : poly_factorize" << endl;
926 #endif
927 
928  poly::get_variables(BHEAD vector<WORD*>(1,argin), with_arghead, true);
929  poly den(BHEAD 0);
930  poly a(poly::argument_to_poly(BHEAD argin, with_arghead, true, &den));
931 
932  // check for modulus calculus
933  WORD modp=poly_determine_modulus(BHEAD true, is_fun_arg, "polynomial factorization");
934  a.setmod(modp,1);
935 
936  // factorize
937  factorized_poly f(polyfact::factorize(a));
938  poly_fix_minus_signs(f);
939 
940  poly num(BHEAD 1);
941  for (int i=0; i<(int)f.factor.size(); i++)
942  if (f.factor[i].is_integer())
943  num = f.factor[i];
944 
945  // determine size
946  int len = with_arghead ? ARGHEAD : 0;
947 
948  if (!num.is_one() || !den.is_one()) {
949  len++;
950  len += MaX(ABS(num[num[1]]), den[den[1]])*2+1;
951  len += with_arghead ? ARGHEAD : 1;
952  }
953 
954  for (int i=0; i<(int)f.factor.size(); i++) {
955  if (!f.factor[i].is_integer()) {
956  len += f.power[i] * f.factor[i].size_of_form_notation();
957  len += f.power[i] * (with_arghead ? ARGHEAD : 1);
958  }
959  }
960 
961  len++;
962 
963  if (argout != NULL) {
964  // check size
965  if (len >= AM.MaxTer) {
966  MLOCK(ErrorMessageLock);
967  MesPrint ("ERROR: factorization doesn't fit in a term");
968  MUNLOCK(ErrorMessageLock);
969  Terminate(-1);
970  }
971  }
972  else {
973  // allocate size
974  argout = (WORD*) Malloc1(len*sizeof(WORD), "poly_factorize");
975  }
976 
977  WORD *old_argout = argout;
978 
979  // constant factor
980  if (!num.is_one() || !den.is_one()) {
981  int n = max(ABS(num[num[1]]), ABS(den[den[1]]));
982 
983  if (with_arghead) {
984  *argout++ = ARGHEAD + 2 + 2*n;
985  for (int i=1; i<ARGHEAD; i++)
986  *argout++ = 0;
987  }
988 
989  *argout++ = 2 + 2*n;
990 
991  for (int i=0; i<n; i++)
992  *argout++ = i<ABS(num[num[1]]) ? num[2+AN.poly_num_vars+i] : 0;
993  for (int i=0; i<n; i++)
994  *argout++ = i<ABS(den[den[1]]) ? den[2+AN.poly_num_vars+i] : 0;
995 
996  *argout++ = SGN(num[num[1]]) * (2*n+1);
997 
998  if (!with_arghead)
999  *argout++ = 0;
1000  else {
1001  if (ToFast(old_argout, old_argout))
1002  argout = old_argout+2;
1003  }
1004  }
1005 
1006  // non-constant factors
1007  for (int i=0; i<(int)f.factor.size(); i++)
1008  if (!f.factor[i].is_integer())
1009  for (int j=0; j<f.power[i]; j++) {
1010  poly::poly_to_argument(f.factor[i],argout,with_arghead);
1011 
1012  if (with_arghead)
1013  argout += *argout > 0 ? *argout : 2;
1014  else {
1015  while (*argout!=0) argout+=*argout;
1016  argout++;
1017  }
1018  }
1019 
1020  *argout=0;
1021 
1022  poly_free_poly_vars(BHEAD "AN.poly_vars_factorize");
1023 
1024  // reset modulo calculation
1025  AN.ncmod = AC.ncmod;
1026 
1027  return old_argout;
1028 }
1029 
1030 /*
1031  #] poly_factorize :
1032  #[ poly_factorize_argument :
1033 */
1034 
1047 int poly_factorize_argument(PHEAD WORD *argin, WORD *argout) {
1048 
1049 #ifdef DEBUG
1050  cout << "*** [" << thetime() << "] CALL : poly_factorize_argument" << endl;
1051 #endif
1052 
1053  poly_factorize(BHEAD argin,argout,true,true);
1054  return 0;
1055 }
1056 
1057 /*
1058  #] poly_factorize_argument :
1059  #[ poly_factorize_dollar :
1060 */
1061 
1074 WORD *poly_factorize_dollar (PHEAD WORD *argin) {
1075 
1076 #ifdef DEBUG
1077  cout << "*** [" << thetime() << "] CALL : poly_factorize_dollar" << endl;
1078 #endif
1079 
1080  return poly_factorize(BHEAD argin,NULL,false,false);
1081 }
1082 
1083 /*
1084  #] poly_factorize_dollar :
1085  #[ poly_factorize_expression :
1086 */
1087 
1101 
1102 #ifdef DEBUG
1103  cout << "*** [" << thetime() << "] CALL : poly_factorize_expression" << endl;
1104 #endif
1105 
1106  GETIDENTITY;
1107 
1108  if (AT.WorkPointer + AM.MaxTer > AT.WorkTop ) {
1109  MLOCK(ErrorMessageLock);
1110  MesWork();
1111  MUNLOCK(ErrorMessageLock);
1112  Terminate(-1);
1113  }
1114 
1115  WORD *term = AT.WorkPointer;
1116  WORD startebuf = cbuf[AT.ebufnum].numrhs;
1117  FILEHANDLE *file;
1118  POSITION pos;
1119 
1120  FILEHANDLE *oldinfile = AR.infile;
1121  FILEHANDLE *oldoutfile = AR.outfile;
1122  WORD oldBracketOn = AR.BracketOn;
1123  WORD *oldBrackBuf = AT.BrackBuf;
1124  WORD oldbracketindexflag = AT.bracketindexflag;
1125  char oldCommercial[COMMERCIALSIZE+2];
1126 
1127  strcpy(oldCommercial, (char*)AC.Commercial);
1128  strcpy((char*)AC.Commercial, "factorize");
1129 
1130  // locate is the input
1131  if (expr->status == HIDDENGEXPRESSION || expr->status == HIDDENLEXPRESSION ||
1132  expr->status == INTOHIDEGEXPRESSION || expr->status == INTOHIDELEXPRESSION) {
1133  AR.InHiBuf = 0; file = AR.hidefile; AR.GetFile = 2;
1134  }
1135  else {
1136  AR.InInBuf = 0; file = AR.outfile; AR.GetFile = 0;
1137  }
1138 
1139  // read and write to expression file
1140  AR.infile = AR.outfile = file;
1141 
1142  // dummy indices are not allowed
1143  if (expr->numdummies > 0) {
1144  MesPrint("ERROR: factorization with dummy indices not implemented");
1145  Terminate(-1);
1146  }
1147 
1148  // determine whether the expression in on file or in memory
1149  if (file->handle >= 0) {
1150  pos = expr->onfile;
1151  SeekFile(file->handle,&pos,SEEK_SET);
1152  if (ISNOTEQUALPOS(pos,expr->onfile)) {
1153  MesPrint("ERROR: something wrong in scratch file [poly_factorize_expression]");
1154  Terminate(-1);
1155  }
1156  file->POposition = expr->onfile;
1157  file->POfull = file->PObuffer;
1158  if (expr->status == HIDDENGEXPRESSION)
1159  AR.InHiBuf = 0;
1160  else
1161  AR.InInBuf = 0;
1162  }
1163  else {
1164  file->POfill = (WORD *)((UBYTE *)(file->PObuffer)+BASEPOSITION(expr->onfile));
1165  }
1166 
1167  SetScratch(AR.infile, &(expr->onfile));
1168 
1169  // read the first header term
1170  WORD size = GetTerm(BHEAD term);
1171  if (size <= 0) {
1172  MesPrint ("ERROR: something wrong with expression [poly_factorize_expression]");
1173  Terminate(-1);
1174  }
1175 
1176  // store position: this is where the output will go
1177  pos = expr->onfile;
1178  ADDPOS(pos, size*sizeof(WORD));
1179 
1180  // use polynomial as buffer, because it is easy to extend
1181  poly buffer(BHEAD 0);
1182  int bufpos = 0;
1183  int sumcommu = 0;
1184 
1185  // read all terms
1186  while (GetTerm(BHEAD term)) {
1187  // substitute non-symbols by extra symbols
1188  sumcommu += DoesCommu(term);
1189  if ( sumcommu > 1 ) {
1190  MesPrint("ERROR: Cannot factorize an expression with more than one noncommuting object");
1191  Terminate(-1);
1192  }
1193  buffer.check_memory(bufpos);
1194  if (LocalConvertToPoly(BHEAD term, buffer.terms + bufpos, startebuf,0) < 0) {
1195  MesPrint("ERROR: in LocalConvertToPoly [factorize_expression]");
1196  Terminate(-1);
1197  }
1198  bufpos += *(buffer.terms + bufpos);
1199  }
1200  buffer[bufpos] = 0;
1201 
1202  // parse the polynomial
1203 
1204  AN.poly_num_vars = 0;
1205  poly::get_variables(BHEAD vector<WORD*>(1,buffer.terms), false, true);
1206  poly den(BHEAD 0);
1207  poly a(poly::argument_to_poly(BHEAD buffer.terms, false, true, &den));
1208 
1209  // check for modulus calculus
1210  WORD modp=poly_determine_modulus(BHEAD true, false, "polynomial factorization");
1211  a.setmod(modp,1);
1212 
1213  // create output
1214  SetScratch(file, &pos);
1215  NewSort(BHEAD0);
1216 
1217  CBUF *C = cbuf+AC.cbufnum;
1218  CBUF *CC = cbuf+AT.ebufnum;
1219 
1220  // turn brackets on. We force the existence of a bracket index.
1221  WORD nexpr = expr - Expressions;
1222  AR.BracketOn = 1;
1223  AT.BrackBuf = AM.BracketFactors;
1224  AT.bracketindexflag = 1;
1225  ClearBracketIndex(-nexpr-2); // Clears the index made during primary generation
1226  OpenBracketIndex(nexpr); // Set up a new index
1227 
1228  if (a.is_zero()) {
1229  expr->numfactors = 1;
1230  }
1231  else if (a.is_one() && den.is_one()) {
1232  expr->numfactors = 1;
1233 
1234  term[0] = 8;
1235  term[1] = SYMBOL;
1236  term[2] = 4;
1237  term[3] = FACTORSYMBOL;
1238  term[4] = 1;
1239  term[5] = 1;
1240  term[6] = 1;
1241  term[7] = 3;
1242 
1243  AT.WorkPointer += *term;
1244  Generator(BHEAD term, C->numlhs);
1245  AT.WorkPointer = term;
1246  }
1247  else {
1248  factorized_poly fac;
1249  bool iszero = false;
1250 
1251  if (!(expr->vflags & ISFACTORIZED)) {
1252  // factorize the polynomial
1253  fac = polyfact::factorize(a);
1254  poly_fix_minus_signs(fac);
1255  }
1256  else {
1257  int factorsymbol=-1;
1258  for (int i=0; i<AN.poly_num_vars; i++)
1259  if (AN.poly_vars[i] == FACTORSYMBOL)
1260  factorsymbol = i;
1261 
1262  poly denpow(BHEAD 1);
1263 
1264  // already factorized, so factorize the factors
1265  for (int i=1; i<=expr->numfactors; i++) {
1266  poly origfac(a.coefficient(factorsymbol, i));
1267  factorized_poly fac2;
1268  if (origfac.is_zero())
1269  iszero=true;
1270  else {
1271  fac2 = polyfact::factorize(origfac);
1272  poly_fix_minus_signs(fac2);
1273  denpow *= den;
1274  }
1275  for (int j=0; j<(int)fac2.power.size(); j++)
1276  fac.add_factor(fac2.factor[j], fac2.power[j]);
1277  }
1278 
1279  // update denominator, since each factor was scaled
1280  den=denpow;
1281  }
1282 
1283  expr->numfactors = 0;
1284 
1285  // coefficient
1286  poly num(BHEAD 1);
1287  for (int i=0; i<(int)fac.factor.size(); i++)
1288  if (fac.factor[i].is_integer())
1289  num *= fac.factor[i];
1290 
1291  poly gcd(polygcd::integer_gcd(num,den));
1292  den/=gcd;
1293  num/=gcd;
1294 
1295  if (iszero)
1296  expr->numfactors++;
1297 
1298  if (!iszero || (expr->vflags & KEEPZERO)) {
1299  if (!num.is_one() || !den.is_one()) {
1300  expr->numfactors++;
1301 
1302  int n = max(ABS(num[num[1]]), ABS(den[den[1]]));
1303 
1304  term[0] = 6 + 2*n;
1305  term[1] = SYMBOL;
1306  term[2] = 4;
1307  term[3] = FACTORSYMBOL;
1308  term[4] = expr->numfactors;
1309  for (int i=0; i<n; i++) {
1310  term[5+i] = i<ABS(num[num[1]]) ? num[2+AN.poly_num_vars+i] : 0;
1311  term[5+n+i] = i<ABS(den[den[1]]) ? den[2+AN.poly_num_vars+i] : 0;
1312  }
1313  term[5+2*n] = SGN(num[num[1]]) * (2*n+1);
1314  AT.WorkPointer += *term;
1315  Generator(BHEAD term, C->numlhs);
1316  AT.WorkPointer = term;
1317  }
1318 
1319  vector<poly> fac_arg(fac.factor.size(), poly(BHEAD 0));
1320 
1321  // convert the non-constant factors to Form-style arguments
1322  for (int i=0; i<(int)fac.factor.size(); i++)
1323  if (!fac.factor[i].is_integer()) {
1324  buffer.check_memory(fac.factor[i].size_of_form_notation()+1);
1325  poly::poly_to_argument(fac.factor[i], buffer.terms, false);
1326  NewSort(BHEAD0);
1327 
1328  for (WORD *t=buffer.terms; *t!=0; t+=*t) {
1329  // substitute extra symbols
1330  if (ConvertFromPoly(BHEAD t, term, numxsymbol, CC->numrhs-startebuf+numxsymbol,
1331  startebuf-numxsymbol, 1) <= 0 ) {
1332  MesPrint("ERROR: in ConvertFromPoly [factorize_expression]");
1333  Terminate(-1);
1334  return(-1);
1335  }
1336 
1337  // store term
1338  AT.WorkPointer += *term;
1339  Generator(BHEAD term, C->numlhs);
1340  AT.WorkPointer = term;
1341  }
1342 
1343  // sort and store in buffer
1344  WORD *buffer;
1345  if (EndSort(BHEAD (WORD *)((VOID *)(&buffer)),2) < 0) return -1;
1346 
1347  LONG bufsize=0;
1348  for (WORD *t=buffer; *t!=0; t+=*t)
1349  bufsize+=*t;
1350 
1351  fac_arg[i].check_memory(bufsize+ARGHEAD+1);
1352 
1353  for (int j=0; j<ARGHEAD; j++)
1354  fac_arg[i].terms[j] = 0;
1355  fac_arg[i].terms[0] = ARGHEAD + bufsize;
1356  WCOPY(fac_arg[i].terms+ARGHEAD, buffer, bufsize);
1357  M_free(buffer, "polynomial factorization");
1358  }
1359 
1360  // compare and sort the factors in Form notation
1361  vector<int> order;
1362  vector<vector<int> > comp(fac.factor.size(), vector<int>(fac.factor.size(), 0));
1363 
1364  for (int i=0; i<(int)fac.factor.size(); i++)
1365  if (!fac.factor[i].is_integer()) {
1366  order.push_back(i);
1367 
1368  for (int j=i+1; j<(int)fac.factor.size(); j++)
1369  if (!fac.factor[j].is_integer()) {
1370  comp[i][j] = CompArg(fac_arg[j].terms, fac_arg[i].terms);
1371  comp[j][i] = -comp[i][j];
1372  }
1373  }
1374 
1375  for (int i=0; i<(int)order.size(); i++)
1376  for (int j=0; j+1<(int)order.size(); j++)
1377  if (comp[order[i]][order[j]] == 1)
1378  swap(order[i],order[j]);
1379 
1380  // create the final expression
1381  for (int i=0; i<(int)order.size(); i++)
1382  for (int j=0; j<fac.power[order[i]]; j++) {
1383 
1384  expr->numfactors++;
1385 
1386  WORD *tstop = fac_arg[order[i]].terms + *fac_arg[order[i]].terms;
1387  for (WORD *t=fac_arg[order[i]].terms+ARGHEAD; t<tstop; t+=*t) {
1388 
1389  WCOPY(term+4, t, *t);
1390 
1391  // add special symbol "factor_"
1392  *term = *(term+4) + 4;
1393  *(term+1) = SYMBOL;
1394  *(term+2) = 4;
1395  *(term+3) = FACTORSYMBOL;
1396  *(term+4) = expr->numfactors;
1397 
1398  // store term
1399  AT.WorkPointer += *term;
1400  Generator(BHEAD term, C->numlhs);
1401  AT.WorkPointer = term;
1402  }
1403  }
1404  }
1405  }
1406 
1407  // final sorting
1408  if (EndSort(BHEAD NULL,0) < 0) {
1409  LowerSortLevel();
1410  Terminate(-1);
1411  }
1412 
1413  // set factorized flag
1414  if (expr->numfactors > 0)
1415  expr->vflags |= ISFACTORIZED;
1416 
1417  // clean up
1418  AR.infile = oldinfile;
1419  AR.outfile = oldoutfile;
1420  AR.BracketOn = oldBracketOn;
1421  AT.BrackBuf = oldBrackBuf;
1422  AT.bracketindexflag = oldbracketindexflag;
1423  strcpy((char*)AC.Commercial, oldCommercial);
1424 
1425  poly_free_poly_vars(BHEAD "AN.poly_vars_factorize_expression");
1426 
1427  return 0;
1428 }
1429 
1430 /*
1431  #] poly_factorize_expression :
1432  #[ poly_unfactorize_expression :
1433 */
1434 
1447 #if ( SUBEXPSIZE == 5 )
1448 static WORD genericterm[] = {38,1,4,FACTORSYMBOL,0
1449  ,EXPRESSION,15,0,1,0,13,10,8,1,4,FACTORSYMBOL,0,1,1,3
1450  ,EXPRESSION,15,0,1,0,13,10,8,1,4,FACTORSYMBOL,0,1,1,3
1451  ,1,1,3,0};
1452 static WORD genericterm2[] = {23,1,4,FACTORSYMBOL,0
1453  ,EXPRESSION,15,0,1,0,13,10,8,1,4,FACTORSYMBOL,0,1,1,3
1454  ,1,1,3,0};
1455 #endif
1456 
1457 int poly_unfactorize_expression(EXPRESSIONS expr)
1458 {
1459  GETIDENTITY;
1460  int i, j, nfac = expr->numfactors, nfacp, nexpr = expr - Expressions;
1461  int expriszero = 0;
1462 
1463  FILEHANDLE *oldinfile = AR.infile;
1464  FILEHANDLE *oldoutfile = AR.outfile;
1465  char oldCommercial[COMMERCIALSIZE+2];
1466 
1467  WORD *oldworkpointer = AT.WorkPointer;
1468  WORD *term = AT.WorkPointer, *t, *w, size;
1469 
1470  FILEHANDLE *file;
1471  POSITION pos, oldpos;
1472 
1473  WORD oldBracketOn = AR.BracketOn;
1474  WORD *oldBrackBuf = AT.BrackBuf;
1475  CBUF *C = cbuf+AC.cbufnum;
1476 
1477  if ( ( expr->vflags & ISFACTORIZED ) == 0 ) return(0);
1478 
1479  if ( AT.WorkPointer + AM.MaxTer > AT.WorkTop ) {
1480  MLOCK(ErrorMessageLock);
1481  MesWork();
1482  MUNLOCK(ErrorMessageLock);
1483  Terminate(-1);
1484  }
1485 
1486  oldpos = AS.OldOnFile[nexpr];
1487  AS.OldOnFile[nexpr] = expr->onfile;
1488 
1489  strcpy(oldCommercial, (char*)AC.Commercial);
1490  strcpy((char*)AC.Commercial, "unfactorize");
1491 /*
1492  locate the input
1493 */
1494  if ( expr->status == HIDDENGEXPRESSION || expr->status == HIDDENLEXPRESSION ||
1495  expr->status == INTOHIDEGEXPRESSION || expr->status == INTOHIDELEXPRESSION ) {
1496  AR.InHiBuf = 0; file = AR.hidefile; AR.GetFile = 2;
1497  }
1498  else {
1499  AR.InInBuf = 0; file = AR.outfile; AR.GetFile = 0;
1500  }
1501 /*
1502  read and write to expression file
1503 */
1504  AR.infile = AR.outfile = file;
1505 /*
1506  set the input file to the correct position
1507 */
1508  if ( file->handle >= 0 ) {
1509  pos = expr->onfile;
1510  SeekFile(file->handle,&pos,SEEK_SET);
1511  if (ISNOTEQUALPOS(pos,expr->onfile)) {
1512  MesPrint("ERROR: something wrong in scratch file unfactorize_expression");
1513  Terminate(-1);
1514  }
1515  file->POposition = expr->onfile;
1516  file->POfull = file->PObuffer;
1517  if ( expr->status == HIDDENGEXPRESSION )
1518  AR.InHiBuf = 0;
1519  else
1520  AR.InInBuf = 0;
1521  }
1522  else {
1523  file->POfill = (WORD *)((UBYTE *)(file->PObuffer)+BASEPOSITION(expr->onfile));
1524  }
1525  SetScratch(AR.infile, &(expr->onfile));
1526 /*
1527  Test for whether the first factor is zero.
1528 */
1529  if ( GetFirstBracket(term,nexpr) < 0 ) Terminate(-1);
1530  if ( term[4] != 1 || *term != 8 || term[1] != SYMBOL || term[3] != FACTORSYMBOL || term[4] != 1 ) {
1531  expriszero = 1;
1532  }
1533  SetScratch(AR.infile, &(expr->onfile));
1534 /*
1535  Read the prototype. After this we have the file ready for the output at pos.
1536 */
1537  size = GetTerm(BHEAD term);
1538  if ( size <= 0 ) {
1539  MesPrint ("ERROR: something wrong with expression unfactorize_expression");
1540  Terminate(-1);
1541  }
1542  pos = expr->onfile;
1543  ADDPOS(pos, size*sizeof(WORD));
1544 /*
1545  Set the brackets straight
1546 */
1547  AR.BracketOn = 1;
1548  AT.BrackBuf = AM.BracketFactors;
1549  AT.bracketinfo = 0;
1550  while ( nfac > 2 ) {
1551  nfacp = nfac - nfac%2;
1552 /*
1553  Prepare the bracket index. We have:
1554  e->bracketinfo: the old input bracket index
1555  e->newbracketinfo: the bracket index made for our current input
1556  We need to keep e->bracketinfo in case other workers need it (InParallel)
1557  Hence we work with AT.bracketinfo which takes priority.
1558  Note that in Processor we forced a newbracketinfo to be made.
1559 */
1560  if ( AT.bracketinfo != 0 ) ClearBracketIndex(-1);
1561  AT.bracketinfo = expr->newbracketinfo;
1562  OpenBracketIndex(nexpr);
1563 /*
1564  Now emulate the terms:
1565  sum_(i,0,nfacp,2,factor_^(i/2+1)*F[factor_^(i+1)]*F[factor_^(i+2)])
1566  +factor_^(nfacp/2+1)*F[factor_^nfac]
1567 */
1568  NewSort(BHEAD0);
1569  if ( expriszero == 0 ) {
1570  for ( i = 0; i < nfacp; i += 2 ) {
1571  t = genericterm; w = term = oldworkpointer;
1572  j = *t; NCOPY(w,t,j);
1573  term[4] = i/2+1;
1574  term[7] = nexpr;
1575  term[16] = i+1;
1576  term[22] = nexpr;
1577  term[31] = i+2;
1578  AT.WorkPointer = term + *term;
1579  Generator(BHEAD term, C->numlhs);
1580  }
1581  if ( nfac > nfacp ) {
1582  t = genericterm2; w = term = oldworkpointer;
1583  j = *t; NCOPY(w,t,j);
1584  term[4] = i/2+1;
1585  term[7] = nexpr;
1586  term[16] = nfac;
1587  AT.WorkPointer = term + *term;
1588  Generator(BHEAD term, C->numlhs);
1589  }
1590  }
1591  if ( EndSort(BHEAD AM.S0->sBuffer,0) < 0 ) {
1592  LowerSortLevel();
1593  Terminate(-1);
1594  }
1595 /*
1596  Set the file back into reading position
1597 */
1598  SetScratch(file, &pos);
1599  nfac = (nfac+1)/2;
1600  if ( expriszero ) { nfac = 1; }
1601  }
1602  if ( AT.bracketinfo != 0 ) ClearBracketIndex(-1);
1603  AT.bracketinfo = expr->newbracketinfo;
1604  expr->newbracketinfo = 0;
1605 /*
1606  Reset the brackets to make them ready for the final pass
1607 */
1608  AR.BracketOn = oldBracketOn;
1609  AT.BrackBuf = oldBrackBuf;
1610  if ( AR.BracketOn ) OpenBracketIndex(nexpr);
1611 /*
1612  We distinguish two cases: nfac == 2 and nfac == 1
1613  After preparing the term we skip the factor_ part.
1614 */
1615  NewSort(BHEAD0);
1616  if ( expriszero == 0 ) {
1617  if ( nfac == 1 ) {
1618  t = genericterm2; w = term = oldworkpointer;
1619  j = *t; NCOPY(w,t,j);
1620  term[7] = nexpr;
1621  term[16] = nfac;
1622  }
1623  else if ( nfac == 2 ) {
1624  t = genericterm; w = term = oldworkpointer;
1625  j = *t; NCOPY(w,t,j);
1626  term[7] = nexpr;
1627  term[16] = 1;
1628  term[22] = nexpr;
1629  term[31] = 2;
1630  }
1631  else {
1632  AS.OldOnFile[nexpr] = oldpos;
1633  return(-1);
1634  }
1635  term[4] = term[0]-4;
1636  term += 4;
1637  AT.WorkPointer = term + *term;
1638  Generator(BHEAD term, C->numlhs);
1639  }
1640  if ( EndSort(BHEAD AM.S0->sBuffer,0) < 0 ) {
1641  LowerSortLevel();
1642  Terminate(-1);
1643  }
1644 /*
1645  Final Cleanup
1646 */
1647  expr->numfactors = 0;
1648  expr->vflags &= ~ISFACTORIZED;
1649  if ( AT.bracketinfo != 0 ) ClearBracketIndex(-1);
1650 
1651  AR.infile = oldinfile;
1652  AR.outfile = oldoutfile;
1653  strcpy((char*)AC.Commercial, oldCommercial);
1654  AT.WorkPointer = oldworkpointer;
1655  AS.OldOnFile[nexpr] = oldpos;
1656 
1657  return(0);
1658 }
1659 
1660 /*
1661  #] poly_unfactorize_expression :
1662  #[ poly_inverse :
1663 */
1664 
1665 WORD *poly_inverse(PHEAD WORD *arga, WORD *argb) {
1666 
1667 #ifdef DEBUG
1668  cout << "*** [" << thetime() << "] CALL : poly_inverse" << endl;
1669 #endif
1670 
1671  // Extract variables
1672  vector<WORD *> e;
1673  e.push_back(arga);
1674  e.push_back(argb);
1675  poly::get_variables(BHEAD e, false, true);
1676 
1677  if (AN.poly_num_vars > 1) {
1678  MLOCK(ErrorMessageLock);
1679  MesPrint ((char*)"ERROR: multivariate polynomial inverse is generally impossible");
1680  MUNLOCK(ErrorMessageLock);
1681  Terminate(-1);
1682  }
1683 
1684  // Convert to polynomials
1685  poly a(poly::argument_to_poly(BHEAD arga, false, true));
1686  poly b(poly::argument_to_poly(BHEAD argb, false, true));
1687 
1688  // Check for modulus calculus
1689  WORD modp=poly_determine_modulus(BHEAD true, true, "polynomial inverse");
1690  a.setmod(modp,1);
1691  b.setmod(modp,1);
1692 
1693  if (modp == 0) {
1694  vector<int> x(1,0);
1695  modp = polyfact::choose_prime(a.integer_lcoeff()*b.integer_lcoeff(), x);
1696  }
1697 
1698  poly amodp(a,modp,1);
1699  poly bmodp(b,modp,1);
1700 
1701  // Calculate gcd
1702  vector<poly> xgcd(polyfact::extended_gcd_Euclidean_lifted(amodp,bmodp));
1703  poly invamodp(xgcd[0]);
1704  poly invbmodp(xgcd[1]);
1705 
1706  if (!((invamodp * amodp) % bmodp).is_one()) {
1707  MLOCK(ErrorMessageLock);
1708  MesPrint ((char*)"ERROR: polynomial inverse does not exist");
1709  MUNLOCK(ErrorMessageLock);
1710  Terminate(-1);
1711  }
1712 
1713  // estimate of the size of the Form notation; might be extended later
1714  int ressize = invamodp.size_of_form_notation()+1;
1715  WORD *res = (WORD *)Malloc1(ressize*sizeof(WORD), "poly_inverse");
1716 
1717  // initialize polynomials to store the result
1718  poly primepower(BHEAD modp);
1719  poly inva(invamodp,modp,1);
1720  poly invb(invbmodp,modp,1);
1721 
1722  while (true) {
1723  // convert to Form notation
1724  int j=0;
1725  WORD n=0;
1726  for (int i=1; i<inva[0]; i+=inva[i]) {
1727 
1728  // check whether res should be extended
1729  while (ressize < j + 2*ABS(inva[i+inva[i]-1]) + (inva[i+1]>0?4:0) + 3) {
1730  int newressize = 2*ressize;
1731 
1732  WORD *newres = (WORD *)Malloc1(newressize*sizeof(WORD), "poly_inverse");
1733  WCOPY(newres, res, ressize);
1734  M_free(res, "poly_inverse");
1735  res = newres;
1736  ressize = newressize;
1737  }
1738 
1739  res[j] = 1;
1740  if (inva[i+1]>0) {
1741  res[j+res[j]++] = SYMBOL;
1742  res[j+res[j]++] = 4;
1743  res[j+res[j]++] = AN.poly_vars[0];
1744  res[j+res[j]++] = inva[i+1];
1745  }
1746  MakeLongRational(BHEAD (UWORD *)&inva[i+2], inva[i+inva[i]-1],
1747  (UWORD*)&primepower.terms[3], primepower.terms[primepower.terms[1]],
1748  (UWORD *)&res[j+res[j]], &n);
1749  res[j] += 2*ABS(n);
1750  res[j+res[j]++] = SGN(n)*(2*ABS(n)+1);
1751  j += res[j];
1752  }
1753  res[j]=0;
1754 
1755  // if modulus calculus is set, this is the answer
1756  if (a.modp != 0) break;
1757 
1758  // otherwise check over integers
1759  poly den(BHEAD 0);
1760  poly check(poly::argument_to_poly(BHEAD res, false, true, &den));
1761  if (poly::divides(b.integer_lcoeff(), check.integer_lcoeff())) {
1762  check = check*a - den;
1763  if (poly::divides(b, check)) break;
1764  }
1765 
1766  // if incorrect, lift with quadratic p-adic Newton's iteration.
1767  poly error((poly(BHEAD 1) - a*inva - b*invb) / primepower);
1768  poly errormodpp(error, modp, inva.modn);
1769 
1770  inva.modn *= 2;
1771  invb.modn *= 2;
1772 
1773  poly dinva((inva * errormodpp) % b);
1774  poly dinvb((invb * errormodpp) % a);
1775 
1776  inva += dinva * primepower;
1777  invb += dinvb * primepower;
1778 
1779  primepower *= primepower;
1780  }
1781 
1782  // clean up and reset modulo calculation
1783  poly_free_poly_vars(BHEAD "AN.poly_vars_inverse");
1784 
1785  AN.ncmod = AC.ncmod;
1786  return res;
1787 }
1788 
1789 /*
1790  #] poly_inverse :
1791  #[ poly_mul :
1792 */
1793 
1794 WORD *poly_mul(PHEAD WORD *a, WORD *b) {
1795 
1796 #ifdef DEBUG
1797  cout << "*** [" << thetime() << "] CALL : poly_mul" << endl;
1798 #endif
1799 
1800  // Extract variables
1801  vector<WORD *> e;
1802  e.push_back(a);
1803  e.push_back(b);
1804  poly::get_variables(BHEAD e, false, false); // TODO: any performance effect by sort_vars=true?
1805 
1806  // Convert to polynomials
1807  poly dena(BHEAD 0);
1808  poly denb(BHEAD 0);
1809  poly pa(poly::argument_to_poly(BHEAD a, false, true, &dena));
1810  poly pb(poly::argument_to_poly(BHEAD b, false, true, &denb));
1811 
1812  // Check for modulus calculus
1813  WORD modp = poly_determine_modulus(BHEAD true, true, "polynomial multiplication");
1814  pa.setmod(modp, 1);
1815 
1816  // NOTE: mul_ is currently implemented by translating negative powers of
1817  // symbols to extra symbols. For future improvement, it may be better to
1818  // compute
1819  // (content(a) * content(b)) * (a/content(a)) * (b/content(b))
1820  // to avoid introducing extra symbols for "mixed" cases, e.g.,
1821  // (1+x) * (1/x) -> (1+x) * (1+Z1_).
1822  assert(dena.is_integer());
1823  assert(denb.is_integer());
1824  assert(modp == 0 || dena.is_one());
1825  assert(modp == 0 || denb.is_one());
1826 
1827  // multiplication
1828  pa *= pb;
1829 
1830  // convert to Form notation
1831  WORD *res;
1832  if (dena.is_one() && denb.is_one()) {
1833  res = (WORD *)Malloc1((pa.size_of_form_notation() + 1) * sizeof(WORD), "poly_mul");
1834  poly::poly_to_argument(pa, res, false);
1835  } else {
1836  dena *= denb;
1837  res = (WORD *)Malloc1((pa.size_of_form_notation_with_den(dena[dena[1]]) + 1) * sizeof(WORD), "poly_mul");
1838  poly::poly_to_argument_with_den(pa, dena[dena[1]], (const UWORD *)&dena[2+AN.poly_num_vars], res, false);
1839  }
1840 
1841  // clean up and reset modulo calculation
1842  poly_free_poly_vars(BHEAD "AN.poly_vars_mul");
1843  AN.ncmod = AC.ncmod;
1844 
1845  return res;
1846 }
1847 
1848 /*
1849  #] poly_mul :
1850  #[ poly_free_poly_vars :
1851 */
1852 
1853 void poly_free_poly_vars(PHEAD const char *text)
1854 {
1855  if ( AN.poly_vars_type == 0 ) {
1856  TermFree(AN.poly_vars, text);
1857  }
1858  else {
1859  M_free(AN.poly_vars, text);
1860  }
1861  AN.poly_num_vars = 0;
1862  AN.poly_vars = 0;
1863 }
1864 
1865 /*
1866  #] poly_free_poly_vars :
1867 */
WORD * poly_gcd(PHEAD WORD *a, WORD *b, WORD fit)
Definition: polywrap.cc:124
WORD Compare1(WORD *, WORD *, WORD)
Definition: sort.c:2536
WORD CompareSymbols(WORD *, WORD *, WORD)
Definition: sort.c:2976
Definition: structs.h:633
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
Definition: notation.c:510
int SymbolNormalize(WORD *)
Definition: normal.c:5014
void poly_sort(PHEAD WORD *a)
Definition: polywrap.cc:557
int poly_ratfun_normalize(PHEAD WORD *term)
Definition: polywrap.cc:719
int poly_factorize_argument(PHEAD WORD *argin, WORD *argout)
Definition: polywrap.cc:1047
void poly_ratfun_read(WORD *a, poly &num, poly &den)
Definition: polywrap.cc:471
Definition: structs.h:938
WORD StoreTerm(PHEAD WORD *)
Definition: sort.c:4333
int poly_factorize_expression(EXPRESSIONS expr)
Definition: polywrap.cc:1100
Definition: poly.h:49
WORD * poly_ratfun_add(PHEAD WORD *t1, WORD *t2)
Definition: polywrap.cc:600
VOID LowerSortLevel()
Definition: sort.c:4727
WORD NewSort(PHEAD0)
Definition: sort.c:592
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:3101
WORD * poly_factorize(PHEAD WORD *argin, WORD *argout, bool with_arghead, bool is_fun_arg)
Definition: polywrap.cc:922
WORD poly_determine_modulus(PHEAD bool multi_error, bool is_fun_arg, string message)
Definition: polywrap.cc:79
static void divmod_heap(const poly &, const poly &, poly &, poly &, bool, bool, bool &)
Definition: poly.cc:1575
WORD * poly_factorize_dollar(PHEAD WORD *argin)
Definition: polywrap.cc:1074
int handle
Definition: structs.h:661
LONG EndSort(PHEAD WORD *, int)
Definition: sort.c:682