52const int POLYWRAP_DENOMPOWER_INCREASE_FACTOR = 2;
81 if (AC.ncmod==0)
return 0;
83 if (!is_fun_arg || (AC.modmode & ALSOFUNARGS)) {
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);
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);
127 cout <<
"*** [" << thetime() <<
"] CALL : poly_gcd" << endl;
151 poly::get_variables(BHEAD e,
false,
true);
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);
161 poly gcd(polygcd::gcd(pa,pb));
164 int newsize = (gcd.size_of_form_notation()+1);
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);
173 res = TermMalloc(
"poly_gcd");
176 res = (WORD *)Malloc1(newsize*
sizeof(WORD),
"poly_gcd");
178 poly::poly_to_argument(gcd, res,
false);
180 poly_free_poly_vars(BHEAD
"AN.poly_vars_qcd");
194WORD *poly_divmod(PHEAD WORD *a, WORD *b,
int divmod, WORD fit) {
197 cout <<
"*** [" << thetime() <<
"] CALL : poly_divmod" << endl;
207 poly::get_variables(BHEAD e,
false,
false);
210 const int DENOMSYMBOL = MAXPOSITIVE;
219 AN.poly_vars[AN.poly_num_vars++] = DENOMSYMBOL;
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);
228 poly numres(polygcd::integer_content(pa));
229 poly denres(polygcd::integer_content(pb));
241 poly gcdres(polygcd::integer_gcd(numres,denres));
246 poly lcoeffb(pb.integer_lcoeff());
250 if (!lcoeffb.is_one()) {
252 if (AN.poly_num_vars > 2) {
257 int denompower = 0, prevdenompower = 0;
262 bool div_fail =
true;
265 if(denompower < prevdenompower)
268 MLOCK(ErrorMessageLock);
269 MesPrint ((
char*)
"ERROR: pseudo-division failed in poly_divmod (denompower > INT_MAX)");
270 MUNLOCK(ErrorMessageLock);
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;
292 prevdenompower = denompower;
293 denompower = (denompower==0 ? 1 : denompower*POLYWRAP_DENOMPOWER_INCREASE_FACTOR+1 );
297 pres = (divmod==0 ? pq : pr);
303 int denompower = MaX(0, pa.degree(0) - pb.degree(0) + 1);
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;
315 pres = (divmod==0 ? pa/pb : pa%pb);
322 pres = (divmod==0 ? pa/pb : pa%pb);
334 res = TermMalloc(
"poly_divmod");
337 res = (WORD *)Malloc1(
sizeof(WORD),
"poly_divmod");
345 UWORD *den = (UWORD *)NumberMalloc(
"poly_divmod");
349 int ressize = pres.size_of_form_notation() + pres.number_of_terms()*2*ABS(denres[denres[1]]) + 1;
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);
358 res = TermMalloc(
"poly_divmod");
361 res = (WORD *)Malloc1(ressize*
sizeof(WORD),
"poly_divmod");
365 for (
int i=1; i!=pres[0]; i+=pres[i]) {
370 for (
int j=0; j<AN.poly_num_vars; j++)
371 if (pres[i+1+j] > 0) {
377 res[L+1+res[L+2]++] = AN.poly_vars[j];
378 res[L+1+res[L+2]++] = pres[i+1+j];
381 if (!first) res[L] += res[L+2];
384 WORD nnum = pres[i+pres[i]-1];
385 WCOPY(&res[L+res[L]], &pres[i+pres[i]-1-ABS(nnum)], ABS(nnum));
388 nden = denres[denres[1]];
389 WCOPY(den, &denres[2+AN.poly_num_vars], ABS(nden));
391 if (nden!=1 || den[0]!=1)
392 Simplify(BHEAD (UWORD *)&res[L+res[L]], &nnum, den, &nden);
393 Pack((UWORD *)&res[L+res[L]], &nnum, den, nden);
394 res[L] += 2*ABS(nnum)+1;
395 res[L+res[L]-1] = SGN(nnum)*(2*ABS(nnum)+1);
401 NumberFree(den,
"poly_divmod");
405 poly_free_poly_vars(BHEAD
"AN.poly_vars_divmod");
424WORD *poly_div(PHEAD WORD *a, WORD *b, WORD fit) {
427 cout <<
"*** [" << thetime() <<
"] CALL : poly_div" << endl;
430 return poly_divmod(BHEAD a, b, 0, fit);
445WORD *poly_rem(PHEAD WORD *a, WORD *b, WORD fit) {
448 cout <<
"*** [" << thetime() <<
"] CALL : poly_rem" << endl;
451 return poly_divmod(BHEAD a, b, 1, fit);
474 cout <<
"*** [" << thetime() <<
"] CALL : poly_ratfun_read" << endl;
477 POLY_GETIDENTITY(num);
481 WORD *astop = a+a[1];
483 bool clean = (a[2] & MUSTCLEANPRF) == 0;
487 MLOCK(ErrorMessageLock);
488 MesPrint ((
char*)
"ERROR: PolyRatFun cannot have zero arguments");
489 MUNLOCK(ErrorMessageLock);
493 poly den_num(BHEAD 1),den_den(BHEAD 1);
495 num = poly::argument_to_poly(BHEAD a,
true, !clean, &den_num);
500 den = poly::argument_to_poly(BHEAD a,
true, !clean, &den_den);
505 den =
poly(BHEAD 1, modp, 1);
509 MLOCK(ErrorMessageLock);
510 MesPrint ((
char*)
"ERROR: PolyRatFun cannot have more than two arguments");
511 MUNLOCK(ErrorMessageLock);
516 vector<WORD> minpower(AN.poly_num_vars, MAXPOSITIVE);
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]);
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];
535 poly gcd = polygcd::gcd(num,den);
560 cout <<
"*** [" << thetime() <<
"] CALL : poly_sort" << endl;
562 if (
NewSort(BHEAD0)) { Terminate(-1); }
565 for (
int i=ARGHEAD; i<a[0]; i+=a[i]) {
573 if (
EndSort(BHEAD a+ARGHEAD,1) < 0) {
602 if ( AR.PolyFunExp == 1 )
return PolyRatFunSpecial(BHEAD t1, t2);
605 cout <<
"*** [" << thetime() <<
"] CALL : poly_ratfun_add" << endl;
608 WORD *oldworkpointer = AT.WorkPointer;
613 for (WORD *t=t1+FUNHEAD; t<t1+t1[1];) {
617 for (WORD *t=t2+FUNHEAD; t<t2+t2[1];) {
622 poly::get_variables(BHEAD e,
true,
true);
628 poly num1(BHEAD 0,modp,1), den1(BHEAD 0,modp,1), num2(BHEAD 0,modp,1), den2(BHEAD 0,modp,1);
633 poly num(BHEAD 0),den(BHEAD 0),gcd(BHEAD 0);
637 gcd = polygcd::gcd(den1,den2);
639 num = num1*(den2/gcd) + num2*(den1/gcd);
640 den = (den1/gcd)*den2;
641 gcd = polygcd::gcd(num,den);
644 num = num1*(den2/gcd) + num2*den;
646 gcd = polygcd::gcd(num,gcd);
652 gcd = polygcd::gcd(num,den);
659 if (den.sign() == -1) { num*=
poly(BHEAD -1); den*=
poly(BHEAD -1); }
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);
672 WORD *t = oldworkpointer;
679 poly::poly_to_argument(num,t,
true);
680 if (*t>0 && t[1]==DIRTYFLAG)
682 t += (*t>0 ? *t : 2);
683 poly::poly_to_argument(den,t,
true);
684 if (*t>0 && t[1]==DIRTYFLAG)
686 t += (*t>0 ? *t : 2);
688 oldworkpointer[1] = t - oldworkpointer;
691 poly_free_poly_vars(BHEAD
"AN.poly_vars_ratfun_add");
696 return oldworkpointer;
722 cout <<
"*** [" << thetime() <<
"] CALL : poly_ratfun_normalize" << endl;
726 WORD *tstop = term + *term;
727 int ncoeff = tstop[-1];
728 tstop -= ABS(ncoeff);
731 int num_polyratfun = 0;
733 for (WORD *t=term+1; t<tstop; t+=t[1])
734 if (*t == AR.PolyFun) {
736 if ((t[2] & MUSTCLEANPRF) != 0)
737 num_polyratfun = INT_MAX;
738 if (num_polyratfun > 1)
break;
741 if (num_polyratfun <= 1)
return 0;
743 WORD oldsorttype = AR.SortType;
744 AR.SortType = SORTHIGHFIRST;
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; }
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];) {
766 poly::get_variables(BHEAD e,
true,
true);
773 poly num1(BHEAD (UWORD *)tstop, ncoeff/2, modp, 1);
774 poly den1(BHEAD (UWORD *)tstop+ABS(ncoeff/2), ABS(ncoeff)/2, modp, 1);
778 for (WORD *t=term+1; t<tstop;)
779 if (*t == AR.PolyFun) {
781 poly num2(BHEAD 0,modp,1);
782 poly den2(BHEAD 0,modp,1);
784 if ((t[2] & MUSTCLEANPRF) != 0) {
785 poly gcd1(polygcd::gcd(num2,den2));
790 poly gcd1(polygcd::gcd(num1,den2));
791 poly gcd2(polygcd::gcd(num2,den1));
793 num1 = (num1 / gcd1) * (num2 / gcd2);
794 den1 = (den1 / gcd2) * (den2 / gcd1);
798 if ( s != t ) { NCOPY(s,t,i) }
799 else { t += i; s += i; }
803 if (den1.sign() == -1) { num1*=
poly(BHEAD -1); den1*=
poly(BHEAD -1); }
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);
819 *t++ &= ~MUSTCLEANPRF;
821 poly::poly_to_argument(num1,t,
true);
822 if (*t>0 && t[1]==DIRTYFLAG)
824 t += (*t>0 ? *t : 2);
825 poly::poly_to_argument(den1,t,
true);
826 if (*t>0 && t[1]==DIRTYFLAG)
828 t += (*t>0 ? *t : 2);
838 poly_free_poly_vars(BHEAD
"AN.poly_vars_ratfun_normalize");
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;
848 AR.SortType = oldsorttype;
859 if ( a.factor.empty() )
return;
861 POLY_GETIDENTITY(a.factor[0]);
863 int overall_sign = +1;
866 for (
int i=0; i<(int)a.factor.size(); i++) {
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];
877 if (pow>0 && (var>maxvar || (var==maxvar && pow>maxpow))) {
880 sign = SGN(*(t+*t-1));
886 a.factor[i] *=
poly(BHEAD sign);
887 if (a.power[i] % 2 == 1) overall_sign*=-1;
892 if (overall_sign == -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);
901 a.add_factor(
poly(BHEAD -1), 1);
922WORD *
poly_factorize (PHEAD WORD *argin, WORD *argout,
bool with_arghead,
bool is_fun_arg) {
925 cout <<
"*** [" << thetime() <<
"] CALL : poly_factorize" << endl;
928 poly::get_variables(BHEAD vector<WORD*>(1,argin), with_arghead,
true);
930 poly a(poly::argument_to_poly(BHEAD argin, with_arghead,
true, &den));
938 poly_fix_minus_signs(f);
941 for (
int i=0; i<(int)f.factor.size(); i++)
942 if (f.factor[i].is_integer())
946 int len = with_arghead ? ARGHEAD : 0;
948 if (!num.is_one() || !den.is_one()) {
950 len += MaX(ABS(num[num[1]]), den[den[1]])*2+1;
951 len += with_arghead ? ARGHEAD : 1;
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);
963 if (argout != NULL) {
965 if (len >= AM.MaxTer) {
966 MLOCK(ErrorMessageLock);
967 MesPrint (
"ERROR: factorization doesn't fit in a term");
968 MUNLOCK(ErrorMessageLock);
974 argout = (WORD*) Malloc1(len*
sizeof(WORD),
"poly_factorize");
977 WORD *old_argout = argout;
980 if (!num.is_one() || !den.is_one()) {
981 int n = max(ABS(num[num[1]]), ABS(den[den[1]]));
984 *argout++ = ARGHEAD + 2 + 2*n;
985 for (
int i=1; i<ARGHEAD; i++)
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;
996 *argout++ = SGN(num[num[1]]) * (2*n+1);
1001 if (ToFast(old_argout, old_argout))
1002 argout = old_argout+2;
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);
1013 argout += *argout > 0 ? *argout : 2;
1015 while (*argout!=0) argout+=*argout;
1022 poly_free_poly_vars(BHEAD
"AN.poly_vars_factorize");
1025 AN.ncmod = AC.ncmod;
1050 cout <<
"*** [" << thetime() <<
"] CALL : poly_factorize_argument" << endl;
1077 cout <<
"*** [" << thetime() <<
"] CALL : poly_factorize_dollar" << endl;
1103 cout <<
"*** [" << thetime() <<
"] CALL : poly_factorize_expression" << endl;
1108 if (AT.WorkPointer + AM.MaxTer > AT.WorkTop ) {
1109 MLOCK(ErrorMessageLock);
1111 MUNLOCK(ErrorMessageLock);
1115 WORD *term = AT.WorkPointer;
1116 WORD startebuf = cbuf[AT.ebufnum].numrhs;
1122 WORD oldBracketOn = AR.BracketOn;
1123 WORD *oldBrackBuf = AT.BrackBuf;
1124 WORD oldbracketindexflag = AT.bracketindexflag;
1125 char oldCommercial[COMMERCIALSIZE+2];
1127 strcpy(oldCommercial, (
char*)AC.Commercial);
1128 strcpy((
char*)AC.Commercial,
"factorize");
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;
1136 AR.InInBuf = 0; file = AR.outfile; AR.GetFile = 0;
1140 AR.infile = AR.outfile = file;
1143 if (expr->numdummies > 0) {
1144 MesPrint(
"ERROR: factorization with dummy indices not implemented");
1151 SeekFile(file->
handle,&pos,SEEK_SET);
1152 if (ISNOTEQUALPOS(pos,expr->onfile)) {
1153 MesPrint(
"ERROR: something wrong in scratch file [poly_factorize_expression]");
1156 file->POposition = expr->onfile;
1157 file->POfull = file->PObuffer;
1158 if (expr->status == HIDDENGEXPRESSION)
1164 file->POfill = (WORD *)((UBYTE *)(file->PObuffer)+BASEPOSITION(expr->onfile));
1167 SetScratch(AR.infile, &(expr->onfile));
1170 WORD size = GetTerm(BHEAD term);
1172 MesPrint (
"ERROR: something wrong with expression [poly_factorize_expression]");
1178 ADDPOS(pos, size*
sizeof(WORD));
1181 poly buffer(BHEAD 0);
1186 while (GetTerm(BHEAD term)) {
1188 sumcommu += DoesCommu(term);
1189 if ( sumcommu > 1 ) {
1190 MesPrint(
"ERROR: Cannot factorize an expression with more than one noncommuting object");
1193 buffer.check_memory(bufpos);
1195 MesPrint(
"ERROR: in LocalConvertToPoly [factorize_expression]");
1198 bufpos += *(buffer.terms + bufpos);
1204 AN.poly_num_vars = 0;
1205 poly::get_variables(BHEAD vector<WORD*>(1,buffer.terms),
false,
true);
1207 poly a(poly::argument_to_poly(BHEAD buffer.terms,
false,
true, &den));
1214 SetScratch(file, &pos);
1217 CBUF *C = cbuf+AC.cbufnum;
1218 CBUF *CC = cbuf+AT.ebufnum;
1221 WORD nexpr = expr - Expressions;
1223 AT.BrackBuf = AM.BracketFactors;
1224 AT.bracketindexflag = 1;
1225 ClearBracketIndex(-nexpr-2);
1226 OpenBracketIndex(nexpr);
1229 expr->numfactors = 1;
1231 else if (a.is_one() && den.is_one()) {
1232 expr->numfactors = 1;
1237 term[3] = FACTORSYMBOL;
1243 AT.WorkPointer += *term;
1245 AT.WorkPointer = term;
1249 bool iszero =
false;
1251 if (!(expr->vflags & ISFACTORIZED)) {
1253 fac = polyfact::factorize(a);
1254 poly_fix_minus_signs(fac);
1257 int factorsymbol=-1;
1258 for (
int i=0; i<AN.poly_num_vars; i++)
1259 if (AN.poly_vars[i] == FACTORSYMBOL)
1262 poly denpow(BHEAD 1);
1265 for (
int i=1; i<=expr->numfactors; i++) {
1266 poly origfac(a.coefficient(factorsymbol, i));
1268 if (origfac.is_zero())
1271 fac2 = polyfact::factorize(origfac);
1272 poly_fix_minus_signs(fac2);
1275 for (
int j=0; j<(int)fac2.power.size(); j++)
1276 fac.add_factor(fac2.factor[j], fac2.power[j]);
1283 expr->numfactors = 0;
1287 for (
int i=0; i<(int)fac.factor.size(); i++)
1288 if (fac.factor[i].is_integer())
1289 num *= fac.factor[i];
1291 poly gcd(polygcd::integer_gcd(num,den));
1298 if (!iszero || (expr->vflags & KEEPZERO)) {
1299 if (!num.is_one() || !den.is_one()) {
1302 int n = max(ABS(num[num[1]]), ABS(den[den[1]]));
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;
1313 term[5+2*n] = SGN(num[num[1]]) * (2*n+1);
1314 AT.WorkPointer += *term;
1316 AT.WorkPointer = term;
1319 vector<poly> fac_arg(fac.factor.size(),
poly(BHEAD 0));
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);
1328 for (WORD *t=buffer.terms; *t!=0; t+=*t) {
1330 if (ConvertFromPoly(BHEAD t, term, numxsymbol, CC->numrhs-startebuf+numxsymbol,
1331 startebuf-numxsymbol, 1) <= 0 ) {
1332 MesPrint(
"ERROR: in ConvertFromPoly [factorize_expression]");
1338 AT.WorkPointer += *term;
1340 AT.WorkPointer = term;
1345 if (
EndSort(BHEAD (WORD *)((VOID *)(&buffer)),2) < 0)
return -1;
1348 for (WORD *t=buffer; *t!=0; t+=*t)
1351 fac_arg[i].check_memory(bufsize+ARGHEAD+1);
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");
1362 vector<vector<int> > comp(fac.factor.size(), vector<int>(fac.factor.size(), 0));
1364 for (
int i=0; i<(int)fac.factor.size(); i++)
1365 if (!fac.factor[i].is_integer()) {
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];
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]);
1381 for (
int i=0; i<(int)order.size(); i++)
1382 for (
int j=0; j<fac.power[order[i]]; j++) {
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) {
1389 WCOPY(term+4, t, *t);
1392 *term = *(term+4) + 4;
1395 *(term+3) = FACTORSYMBOL;
1396 *(term+4) = expr->numfactors;
1399 AT.WorkPointer += *term;
1401 AT.WorkPointer = term;
1408 if (
EndSort(BHEAD NULL,0) < 0) {
1414 if (expr->numfactors > 0)
1415 expr->vflags |= ISFACTORIZED;
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);
1425 poly_free_poly_vars(BHEAD
"AN.poly_vars_factorize_expression");
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
1452static WORD genericterm2[] = {23,1,4,FACTORSYMBOL,0
1453 ,EXPRESSION,15,0,1,0,13,10,8,1,4,FACTORSYMBOL,0,1,1,3
1460 int i, j, nfac = expr->numfactors, nfacp, nexpr = expr - Expressions;
1465 char oldCommercial[COMMERCIALSIZE+2];
1467 WORD *oldworkpointer = AT.WorkPointer;
1468 WORD *term = AT.WorkPointer, *t, *w, size;
1473 WORD oldBracketOn = AR.BracketOn;
1474 WORD *oldBrackBuf = AT.BrackBuf;
1475 CBUF *C = cbuf+AC.cbufnum;
1477 if ( ( expr->vflags & ISFACTORIZED ) == 0 )
return(0);
1479 if ( AT.WorkPointer + AM.MaxTer > AT.WorkTop ) {
1480 MLOCK(ErrorMessageLock);
1482 MUNLOCK(ErrorMessageLock);
1486 oldpos = AS.OldOnFile[nexpr];
1487 AS.OldOnFile[nexpr] = expr->onfile;
1489 strcpy(oldCommercial, (
char*)AC.Commercial);
1490 strcpy((
char*)AC.Commercial,
"unfactorize");
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;
1499 AR.InInBuf = 0; file = AR.outfile; AR.GetFile = 0;
1504 AR.infile = AR.outfile = file;
1508 if ( file->
handle >= 0 ) {
1510 SeekFile(file->
handle,&pos,SEEK_SET);
1511 if (ISNOTEQUALPOS(pos,expr->onfile)) {
1512 MesPrint(
"ERROR: something wrong in scratch file unfactorize_expression");
1515 file->POposition = expr->onfile;
1516 file->POfull = file->PObuffer;
1517 if ( expr->status == HIDDENGEXPRESSION )
1523 file->POfill = (WORD *)((UBYTE *)(file->PObuffer)+BASEPOSITION(expr->onfile));
1525 SetScratch(AR.infile, &(expr->onfile));
1529 if ( GetFirstBracket(term,nexpr) < 0 ) Terminate(-1);
1530 if ( term[4] != 1 || *term != 8 || term[1] != SYMBOL || term[3] != FACTORSYMBOL || term[4] != 1 ) {
1533 SetScratch(AR.infile, &(expr->onfile));
1537 size = GetTerm(BHEAD term);
1539 MesPrint (
"ERROR: something wrong with expression unfactorize_expression");
1543 ADDPOS(pos, size*
sizeof(WORD));
1548 AT.BrackBuf = AM.BracketFactors;
1550 while ( nfac > 2 ) {
1551 nfacp = nfac - nfac%2;
1560 if ( AT.bracketinfo != 0 ) ClearBracketIndex(-1);
1561 AT.bracketinfo = expr->newbracketinfo;
1562 OpenBracketIndex(nexpr);
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);
1578 AT.WorkPointer = term + *term;
1581 if ( nfac > nfacp ) {
1582 t = genericterm2; w = term = oldworkpointer;
1583 j = *t; NCOPY(w,t,j);
1587 AT.WorkPointer = term + *term;
1591 if (
EndSort(BHEAD AM.S0->sBuffer,0) < 0 ) {
1598 SetScratch(file, &pos);
1600 if ( expriszero ) { nfac = 1; }
1602 if ( AT.bracketinfo != 0 ) ClearBracketIndex(-1);
1603 AT.bracketinfo = expr->newbracketinfo;
1604 expr->newbracketinfo = 0;
1608 AR.BracketOn = oldBracketOn;
1609 AT.BrackBuf = oldBrackBuf;
1610 if ( AR.BracketOn ) OpenBracketIndex(nexpr);
1616 if ( expriszero == 0 ) {
1618 t = genericterm2; w = term = oldworkpointer;
1619 j = *t; NCOPY(w,t,j);
1623 else if ( nfac == 2 ) {
1624 t = genericterm; w = term = oldworkpointer;
1625 j = *t; NCOPY(w,t,j);
1632 AS.OldOnFile[nexpr] = oldpos;
1635 term[4] = term[0]-4;
1637 AT.WorkPointer = term + *term;
1640 if (
EndSort(BHEAD AM.S0->sBuffer,0) < 0 ) {
1647 expr->numfactors = 0;
1648 expr->vflags &= ~ISFACTORIZED;
1649 if ( AT.bracketinfo != 0 ) ClearBracketIndex(-1);
1651 AR.infile = oldinfile;
1652 AR.outfile = oldoutfile;
1653 strcpy((
char*)AC.Commercial, oldCommercial);
1654 AT.WorkPointer = oldworkpointer;
1655 AS.OldOnFile[nexpr] = oldpos;
1665WORD *poly_inverse(PHEAD WORD *arga, WORD *argb) {
1668 cout <<
"*** [" << thetime() <<
"] CALL : poly_inverse" << endl;
1675 poly::get_variables(BHEAD e,
false,
true);
1677 if (AN.poly_num_vars > 1) {
1678 MLOCK(ErrorMessageLock);
1679 MesPrint ((
char*)
"ERROR: multivariate polynomial inverse is generally impossible");
1680 MUNLOCK(ErrorMessageLock);
1685 poly a(poly::argument_to_poly(BHEAD arga,
false,
true));
1686 poly b(poly::argument_to_poly(BHEAD argb,
false,
true));
1695 modp = polyfact::choose_prime(a.integer_lcoeff()*b.integer_lcoeff(), x);
1698 poly amodp(a,modp,1);
1699 poly bmodp(b,modp,1);
1702 vector<poly> xgcd(polyfact::extended_gcd_Euclidean_lifted(amodp,bmodp));
1703 poly invamodp(xgcd[0]);
1704 poly invbmodp(xgcd[1]);
1706 if (!((invamodp * amodp) % bmodp).is_one()) {
1707 MLOCK(ErrorMessageLock);
1708 MesPrint ((
char*)
"ERROR: polynomial inverse does not exist");
1709 MUNLOCK(ErrorMessageLock);
1714 int ressize = invamodp.size_of_form_notation()+1;
1715 WORD *res = (WORD *)Malloc1(ressize*
sizeof(WORD),
"poly_inverse");
1718 poly primepower(BHEAD modp);
1719 poly inva(invamodp,modp,1);
1720 poly invb(invbmodp,modp,1);
1726 for (
int i=1; i<inva[0]; i+=inva[i]) {
1729 while (ressize < j + 2*ABS(inva[i+inva[i]-1]) + (inva[i+1]>0?4:0) + 3) {
1730 int newressize = 2*ressize;
1732 WORD *newres = (WORD *)Malloc1(newressize*
sizeof(WORD),
"poly_inverse");
1733 WCOPY(newres, res, ressize);
1734 M_free(res,
"poly_inverse");
1736 ressize = newressize;
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];
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);
1750 res[j+res[j]++] = SGN(n)*(2*ABS(n)+1);
1756 if (a.modp != 0)
break;
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;
1767 poly error((
poly(BHEAD 1) - a*inva - b*invb) / primepower);
1768 poly errormodpp(error, modp, inva.modn);
1773 poly dinva((inva * errormodpp) % b);
1774 poly dinvb((invb * errormodpp) % a);
1776 inva += dinva * primepower;
1777 invb += dinvb * primepower;
1779 primepower *= primepower;
1783 poly_free_poly_vars(BHEAD
"AN.poly_vars_inverse");
1785 AN.ncmod = AC.ncmod;
1794WORD *poly_mul(PHEAD WORD *a, WORD *b) {
1797 cout <<
"*** [" << thetime() <<
"] CALL : poly_mul" << endl;
1804 poly::get_variables(BHEAD e,
false,
false);
1809 poly pa(poly::argument_to_poly(BHEAD a,
false,
true, &dena));
1810 poly pb(poly::argument_to_poly(BHEAD b,
false,
true, &denb));
1822 assert(dena.is_integer());
1823 assert(denb.is_integer());
1824 assert(modp == 0 || dena.is_one());
1825 assert(modp == 0 || denb.is_one());
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);
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);
1842 poly_free_poly_vars(BHEAD
"AN.poly_vars_mul");
1843 AN.ncmod = AC.ncmod;
1853void poly_free_poly_vars(PHEAD
const char *text)
1855 if ( AN.poly_vars_type == 0 ) {
1856 TermFree(AN.poly_vars, text);
1859 M_free(AN.poly_vars, text);
1861 AN.poly_num_vars = 0;
static void divmod_heap(const poly &, const poly &, poly &, poly &, bool, bool, bool &)
WORD CompareSymbols(WORD *, WORD *, WORD)
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
LONG EndSort(PHEAD WORD *, int)
WORD Generator(PHEAD WORD *, WORD)
WORD StoreTerm(PHEAD WORD *)
WORD Compare1(WORD *, WORD *, WORD)
int SymbolNormalize(WORD *)
int poly_factorize_expression(EXPRESSIONS expr)
WORD poly_determine_modulus(PHEAD bool multi_error, bool is_fun_arg, string message)
WORD * poly_ratfun_add(PHEAD WORD *t1, WORD *t2)
void poly_ratfun_read(WORD *a, poly &num, poly &den)
void poly_sort(PHEAD WORD *a)
WORD * poly_gcd(PHEAD WORD *a, WORD *b, WORD fit)
int poly_ratfun_normalize(PHEAD WORD *term)
WORD * poly_factorize(PHEAD WORD *argin, WORD *argout, bool with_arghead, bool is_fun_arg)
int poly_unfactorize_expression(EXPRESSIONS expr)
int poly_factorize_argument(PHEAD WORD *argin, WORD *argout)
WORD * poly_factorize_dollar(PHEAD WORD *argin)