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
52const int POLYWRAP_DENOMPOWER_INCREASE_FACTOR = 2;
53
54using namespace std;
55
56/*
57 #[ poly_determine_modulus :
58*/
59
79WORD 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
124WORD *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
194WORD *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*/
424WORD *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*/
445WORD *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
471void 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
557void 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;
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
600WORD *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
719int 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
857void 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
922WORD *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
1047int 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
1074WORD *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) {
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 )
1448static 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};
1452static 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
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 ) {
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 ) {
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
1665WORD *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
1794WORD *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
1853void 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*/
Definition: poly.h:49
static void divmod_heap(const poly &, const poly &, poly &, poly &, bool, bool, bool &)
Definition: poly.cc:1575
WORD CompareSymbols(WORD *, WORD *, WORD)
Definition: sort.c:2976
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
Definition: notation.c:510
WORD NewSort(PHEAD0)
Definition: sort.c:592
LONG EndSort(PHEAD WORD *, int)
Definition: sort.c:682
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:3101
WORD StoreTerm(PHEAD WORD *)
Definition: sort.c:4333
WORD Compare1(WORD *, WORD *, WORD)
Definition: sort.c:2536
int SymbolNormalize(WORD *)
Definition: normal.c:5014
int poly_factorize_expression(EXPRESSIONS expr)
Definition: polywrap.cc:1100
WORD poly_determine_modulus(PHEAD bool multi_error, bool is_fun_arg, string message)
Definition: polywrap.cc:79
WORD * poly_ratfun_add(PHEAD WORD *t1, WORD *t2)
Definition: polywrap.cc:600
void poly_ratfun_read(WORD *a, poly &num, poly &den)
Definition: polywrap.cc:471
void poly_sort(PHEAD WORD *a)
Definition: polywrap.cc:557
WORD * poly_gcd(PHEAD WORD *a, WORD *b, WORD fit)
Definition: polywrap.cc:124
int poly_ratfun_normalize(PHEAD WORD *term)
Definition: polywrap.cc:719
WORD * poly_factorize(PHEAD WORD *argin, WORD *argout, bool with_arghead, bool is_fun_arg)
Definition: polywrap.cc:922
int poly_unfactorize_expression(EXPRESSIONS expr)
Definition: polywrap.cc:1457
int poly_factorize_argument(PHEAD WORD *argin, WORD *argout)
Definition: polywrap.cc:1047
WORD * poly_factorize_dollar(PHEAD WORD *argin)
Definition: polywrap.cc:1074
VOID LowerSortLevel()
Definition: sort.c:4727
Definition: structs.h:938
Definition: structs.h:633
int handle
Definition: structs.h:661