62WORD RatioFind(PHEAD WORD *term, WORD *params)
67 WORD *y1, *y2, n1 = 0, n2 = 0;
81 if ( *m == x1 ) { y1 = m; n1 = m[1]; }
82 else if ( *m == x2 ) { y2 = m; n2 = m[1]; }
85 if ( !y1 || !y2 || ( n1 > 0 && n2 > 0 ) )
return(0);
87 if ( y1 > y2 ) { r = y1; y1 = y2; y2 = r; }
88 *y2 = *m; y2[1] = m[1];
90 *y1 = *m; y1[1] = m[1];
93We have to revise the code
for the second
case.
100 while ( y1 > r ) *--y2 = *--y1;
101 *m++ = SUBEXPRESSION;
107 *term += SUBEXPSIZE-4;
111 *m++ = SUBEXPRESSION;
119 *term += SUBEXPSIZE-6;
120 r = m + 6-SUBEXPSIZE;
121 do { *m++ = *r++; }
while ( r < t );
170WORD RatioGen(PHEAD WORD *term, WORD *params, WORD num, WORD level)
178 WORD ncoef, sign = 0;
179 coef = (UWORD *)AT.WorkPointer;
181 tstops[2] = m = t + *t;
185 if ( *t == SUBEXPRESSION && t[2] == num )
break;
189 tstops[1] = t + t[1];
206 i = n1; n1 = n2; n2 = i;
207 i = x1; x1 = x2; x2 = i;
216 AT.WorkPointer = (WORD *)(coef + 1);
218 for ( i = 0; i <= n2; i++ ) {
219 if ( BinomGen(BHEAD term,level,tstops,x1,x3,n2-n1-i,i,sign&i
220 ,coef,ncoef) )
goto RatioCall;
222 if ( Product(coef,&ncoef,j) )
goto RatioCall;
223 if ( Quotient(coef,&ncoef,i+1) )
goto RatioCall;
225 AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
228 AT.WorkPointer = (WORD *)(coef);
238 AT.WorkPointer = (WORD *)(coef + 1);
240 for ( i = 0; i <= j; i++ ) {
241 if ( BinomGen(BHEAD term,level,tstops,x2,x3,n2-n1-i,i,sign&i
242 ,coef,ncoef) )
goto RatioCall;
244 if ( Product(coef,&ncoef,n1+i) )
goto RatioCall;
245 if ( Quotient(coef,&ncoef,i+1) )
goto RatioCall;
246 AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
251 AT.WorkPointer = (WORD *)(coef + 1);
253 for ( i = 0; i <= j; i++ ) {
254 if ( BinomGen(BHEAD term,level,tstops,x1,x3,i-n1,n2-i,sign&(n2-i)
255 ,coef,ncoef) )
goto RatioCall;
257 if ( Product(coef,&ncoef,n2-i) )
goto RatioCall;
258 if ( Quotient(coef,&ncoef,i+1) )
goto RatioCall;
259 AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
262 AT.WorkPointer = (WORD *)(coef);
277 AT.WorkPointer = (WORD *)(coef + 1);
279 for ( i = 0; i <= j; i++ ) {
280 if ( BinomGen(BHEAD term,level,tstops,x1,x3,i-n1,-n2-i,i&1
281 ,coef,ncoef) )
goto RatioCall;
283 if ( Product(coef,&ncoef,n2+i) )
goto RatioCall;
284 if ( Quotient(coef,&ncoef,i+1) )
goto RatioCall;
285 AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
290 AT.WorkPointer = (WORD *)(coef + 1);
292 for ( i = 0; i <= j; i++ ) {
293 if ( BinomGen(BHEAD term,level,tstops,x2,x3,i-n2,-n1-i,n1&1
294 ,coef,ncoef) )
goto RatioCall;
296 if ( Product(coef,&ncoef,n1+i) )
goto RatioCall;
297 if ( Quotient(coef,&ncoef,i+1) )
goto RatioCall;
298 AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
301 AT.WorkPointer = (WORD *)(coef);
306 MLOCK(ErrorMessageLock);
308 MUNLOCK(ErrorMessageLock);
320WORD BinomGen(PHEAD WORD *term, WORD level, WORD **tstops, WORD x1, WORD x2,
321 WORD pow1, WORD pow2, WORD sign, UWORD *coef, WORD ncoef)
327 termout = AT.WorkPointer;
330 do { *t++ = *r++; }
while ( r < tstops[0] );
333 if ( pow1 == 0 ) t--;
334 else { *t++ = 4; *t++ = x1; *t++ = pow1; }
336 else if ( pow1 == 0 ) {
337 *t++ = 4; *t++ = x2; *t++ = pow2;
340 *t++ = 6; *t++ = x1; *t++ = pow1; *t++ = x2; *t++ = pow2;
343 *t++ = ABS(ncoef) + 3;
345 if ( sign ) *t = -*t;
348 for ( k = 0; k < ncoef; k++ ) *t++ = coef[k];
350 do { *t++ = *r++; }
while ( r < tstops[2] );
351 *termout = WORDDIF(t,termout);
353 if ( AT.WorkPointer > AT.WorkTop ) {
354 MLOCK(ErrorMessageLock);
356 MUNLOCK(ErrorMessageLock);
362 MLOCK(ErrorMessageLock);
364 MUNLOCK(ErrorMessageLock);
367 AT.WorkPointer = termout;
395WORD DoSumF1(PHEAD WORD *term, WORD *params, WORD replac, WORD level)
398 WORD *termout, *t, extractbuff = AT.TMbuff;
399 WORD isum, ival, iinc;
404 if ( ( iinc > 0 && params[4] >= ival )
405 || ( iinc < 0 && params[4] <= ival ) ) {
406 isum = (params[4] - ival)/iinc + 1;
409 termout = AT.WorkPointer;
410 AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
411 if ( AT.WorkPointer > AT.WorkTop ) {
412 MLOCK(ErrorMessageLock);
414 MUNLOCK(ErrorMessageLock);
418 while ( *t != SUBEXPRESSION || t[2] != replac || t[4] != extractbuff )
422 if ( params[2] < 0 ) {
423 while ( *t != INDTOIND || t[2] != -params[2] ) t += t[1];
427 while ( *t > SYMTOSUB || t[2] != params[2] ) t += t[1];
433 while ( C->
Buffer[from] ) {
434 if (
InsertTerm(BHEAD term,replac,extractbuff,C->
Buffer+from,termout,0) < 0 )
goto SumF1Call;
435 AT.WorkPointer = termout + *termout;
436 if (
Generator(BHEAD termout,level) < 0 )
goto SumF1Call;
440 }
while ( --isum > 0 );
441 AT.WorkPointer = termout;
444 MLOCK(ErrorMessageLock);
446 MUNLOCK(ErrorMessageLock);
461WORD Glue(PHEAD WORD *term1, WORD *term2, WORD *sub, WORD insert)
465 WORD ncoef, *t, *t1, *t2, i, nc2, nc3, old, newer;
466 coef = (UWORD *)(TermMalloc(
"Glue"));
471 old = WORDDIF(t,term1);
477 while ( --i >= 0 ) *t2++ = *t1++;
479 nc2 = WildFill(BHEAD t,term2,sub);
484 newer = WORDDIF(t,term1);
485 if ( MulRat(BHEAD (UWORD *)t,REDLENG(nc2),coef,ncoef,(UWORD *)t,&nc3) ) {
486 MLOCK(ErrorMessageLock);
488 MUNLOCK(ErrorMessageLock);
489 TermFree(coef,
"Glue");
494 *t++ = (nc3 >= 0)?i:-i;
495 *term1 = WORDDIF(t,term1);
511 TermFree(coef,
"Glue");
520WORD DoSumF2(PHEAD WORD *term, WORD *params, WORD replac, WORD level)
523 WORD *termout, *t, *from, *sub, *to, extractbuff = AT.TMbuff;
524 WORD isum, ival, iinc, insert, i;
528 if ( ( iinc > 0 && params[4] >= ival )
529 || ( iinc < 0 && params[4] <= ival ) ) {
530 isum = (params[4] - ival)/iinc + 1;
533 termout = AT.WorkPointer;
534 AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
535 if ( AT.WorkPointer > AT.WorkTop ) {
536 MLOCK(ErrorMessageLock);
538 MUNLOCK(ErrorMessageLock);
542 while ( *t != SUBEXPRESSION || t[2] != replac || t[4] != extractbuff ) t += t[1];
543 insert = WORDDIF(t,term);
547 while ( from < t ) *to++ = *from++;
550 while ( from < sub ) *to++ = *from++;
556 if ( params[2] < 0 ) {
557 while ( *t != INDTOIND || t[2] != -params[2] ) t += t[1];
561 while ( *t > SYMTOSUB || t[2] != params[2] ) t += t[1];
566 AT.WorkPointer = termout + *termout;
568 if ( ( to + *termout ) > AT.WorkTop ) {
569 MLOCK(ErrorMessageLock);
571 MUNLOCK(ErrorMessageLock);
577 from = AT.WorkPointer;
579 if (
Generator(BHEAD from,level) < 0 )
goto SumF2Call;
580 if ( --isum <= 0 )
break;
583 if ( Glue(BHEAD termout,C->
rhs[replac],sub,insert) < 0 )
goto SumF2Call;
585 AT.WorkPointer = termout;
588 MLOCK(ErrorMessageLock);
590 MUNLOCK(ErrorMessageLock);
609int GCDfunction(PHEAD WORD *term,WORD level)
612 WORD *t, *tstop, *tf, *termout, *tin, *tout, *m, *mnext, *mstop, *mm;
613 int todo, i, ii, j, istart, sign = 1, action = 0;
614 WORD firstshort = 0, firstvalue = 0, gcdisone = 0, mlength, tlength, newlength;
615 WORD totargs = 0, numargs, argsdone = 0, *mh, oldval1, *g, *gcdout = 0;
631 t = term + *term; tlength = t[-1];
632 tstop = t - ABS(tlength);
634 while ( t < tstop ) {
635 if ( *t != GCDFUNCTION ) { t += t[1];
continue; }
636 todo = 1; totargs = 0;
638 while ( tf < t + t[1] ) {
640 if ( *tf > 0 && tf[1] != 0 ) todo = 0;
647 MLOCK(ErrorMessageLock);
648 MesPrint(
"Internal error. Indicated gcd_ function not encountered.");
649 MUNLOCK(ErrorMessageLock);
652 WantAddPointers(totargs);
653 args = AT.pWorkPointer; AT.pWorkPointer += totargs;
667 while ( tf < t + t[1] ) {
668 if ( *tf == -SNUMBER && tf[1] == 0 ) { NEXTARG(tf);
continue; }
669 if ( *tf > 0 || *tf == -DOLLAREXPRESSION || *tf == -EXPRESSION ) {
670 AT.pWorkSpace[args+numargs++] = tf;
671 NEXTARG(tf);
continue;
673 if ( firstshort == 0 ) {
675 if ( *tf <= -FUNCTION ) { firstvalue = -(*tf); }
676 else { firstvalue = tf[1]; }
681 else if ( *tf != firstshort ) {
682 if ( *tf != -INDEX && *tf != -VECTOR && *tf != -MINVECTOR ) {
683 argsdone++; gcdisone = 1;
break;
685 if ( firstshort != -INDEX && firstshort != -VECTOR && firstshort != -MINVECTOR ) {
686 argsdone++; gcdisone = 1;
break;
688 if ( tf[1] != firstvalue ) {
689 argsdone++; gcdisone = 1;
break;
691 if ( *t == -MINVECTOR ) { firstshort = -VECTOR; }
692 if ( firstshort == -MINVECTOR ) { firstshort = -VECTOR; }
694 else if ( *tf > -FUNCTION && *tf != -SNUMBER && tf[1] != firstvalue ) {
695 argsdone++; gcdisone = 1;
break;
697 if ( *tf == -SNUMBER && firstvalue != tf[1] ) {
701 if ( firstvalue == 1 || tf[1] == 1 ) { gcdisone = 1;
break; }
702 if ( firstvalue < 0 && tf[1] < 0 ) {
703 x1 = -firstvalue; x2 = -tf[1]; sign = -1;
706 x1 = ABS(firstvalue); x2 = ABS(tf[1]); sign = 1;
708 while ( ( x3 = x1%x2 ) != 0 ) { x1 = x2; x2 = x3; }
709 firstvalue = ((WORD)x2)*sign;
711 if ( firstvalue == 1 ) { gcdisone = 1;
break; }
715 termout = AT.WorkPointer;
716 AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
717 if ( AT.WorkPointer > AT.WorkTop ) {
718 MLOCK(ErrorMessageLock);
720 MUNLOCK(ErrorMessageLock);
729 i = t - term; tin = term; tout = termout;
731 if ( gcdisone || ( firstshort == -SNUMBER && firstvalue == 1 ) ) {
734 tin += t[1]; tstop = term + *term;
735 while ( tin < tstop ) *tout++ = *tin++;
736 *termout = tout - termout;
737 if ( sign < 0 ) tout[-1] = -tout[-1];
738 AT.WorkPointer = tout;
739 if ( argsdone &&
Generator(BHEAD termout,level) < 0 )
goto CalledFrom;
740 AT.WorkPointer = termout;
741 AT.pWorkPointer = args;
748 if ( numargs == 0 ) {
751 if ( firstshort == 0 )
goto gcdone;
752 if ( firstshort == -SNUMBER ) {
753 *tout++ = SNUMBER; *tout++ = 4; *tout++ = firstvalue; *tout++ = 1;
756 else if ( firstshort == -SYMBOL ) {
757 *tout++ = SYMBOL; *tout++ = 4; *tout++ = firstvalue; *tout++ = 1;
760 else if ( firstshort == -VECTOR || firstshort == -INDEX ) {
761 *tout++ = INDEX; *tout++ = 3; *tout++ = firstvalue;
goto gcdone;
763 else if ( firstshort == -MINVECTOR ) {
765 *tout++ = INDEX; *tout++ = 3; *tout++ = firstvalue;
goto gcdone;
767 else if ( firstshort <= -FUNCTION ) {
768 *tout++ = firstvalue; *tout++ = FUNHEAD; FILLFUN(tout);
772 MLOCK(ErrorMessageLock);
773 MesPrint(
"Internal error. Illegal short argument in GCDfunction.");
774 MUNLOCK(ErrorMessageLock);
786 switch ( firstshort ) {
788 sh[0] = 4; sh[1] = ABS(firstvalue); sh[2] = 1;
789 if ( firstvalue < 0 ) sh[3] = -3;
796 sh[0] = 8; sh[1] = INDEX; sh[2] = 3; sh[3] = firstvalue;
797 sh[4] = 1; sh[5] = 1;
798 if ( firstshort == -MINVECTOR ) sh[6] = -3;
803 sh[0] = 8; sh[1] = SYMBOL; sh[2] = 4; sh[3] = firstvalue; sh[4] = 1;
804 sh[5] = 1; sh[6] = 1; sh[7] = 3; sh[8] = 0;
807 sh[0] = FUNHEAD+4; sh[1] = firstshort; sh[2] = FUNHEAD;
808 for ( i = 2; i < FUNHEAD; i++ ) sh[i+1] = 0;
809 sh[FUNHEAD+1] = 1; sh[FUNHEAD+2] = 1; sh[FUNHEAD+3] = 3; sh[FUNHEAD+4] = 0;
820 for ( i = 1; i < numargs; i++ ) {
821 for ( ii = i; ii > 0; ii-- ) {
822 arg1 = AT.pWorkSpace[args+ii];
823 arg2 = AT.pWorkSpace[args+ii-1];
826 if ( *arg1 == -EXPRESSION )
break;
827 if ( *arg2 == -DOLLAREXPRESSION )
break;
828 AT.pWorkSpace[args+ii] = arg2;
829 AT.pWorkSpace[args+ii-1] = arg1;
833 else if ( *arg2 < 0 ) {
834 AT.pWorkSpace[args+ii] = arg2;
835 AT.pWorkSpace[args+ii-1] = arg1;
838 if ( *arg1 > *arg2 ) {
839 AT.pWorkSpace[args+ii] = arg2;
840 AT.pWorkSpace[args+ii-1] = arg1;
853 for ( i = istart; i < numargs; i++ ) {
854 arg1 = AT.pWorkSpace[args+i];
856 oldval1 = arg1[*arg1]; arg1[*arg1] = 0;
859 GCDterms(BHEAD mh,m,mh); m += *m;
860 if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
861 gcdisone = 1; sign = 1; arg1[*arg1] = oldval1;
goto gcdone;
864 arg1[*arg1] = oldval1;
866 else if ( *arg1 == -DOLLAREXPRESSION ) {
867 if ( ( d = DolToTerms(BHEAD arg1[1]) ) != 0 ) {
870 GCDterms(BHEAD mh,m,mh); m += *m;
872 if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
873 gcdisone = 1; sign = 1;
874 if ( d->factors ) M_free(d->factors,
"Dollar factors");
875 M_free(d,
"Copy of dollar variable");
goto gcdone;
878 if ( d->factors ) M_free(d->factors,
"Dollar factors");
879 M_free(d,
"Copy of dollar variable");
883 mm = CreateExpression(BHEAD arg1[1]);
886 GCDterms(BHEAD mh,m,mh); m += *m;
888 if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
889 gcdisone = 1; sign = 1; M_free(mm,
"CreateExpression");
goto gcdone;
892 M_free(mm,
"CreateExpression");
897 firstshort = -SNUMBER; firstvalue = mh[1] * (mh[3]/3);
899 else if ( mh[1] == SYMBOL ) {
900 firstshort = -SYMBOL; firstvalue = mh[3];
902 else if ( mh[1] == INDEX ) {
903 firstshort = -INDEX; firstvalue = mh[3];
904 if ( mh[6] == -3 ) firstshort = -MINVECTOR;
906 else if ( mh[1] >= FUNCTION ) {
907 firstshort = -mh[1]; firstvalue = mh[1];
928 for ( i = 0; i < numargs; i++ ) {
929 arg1 = AT.pWorkSpace[args+i];
930 if ( *arg1 > 0 && arg1[ARGHEAD]+ARGHEAD == *arg1 ) {
935 arg2 = AT.pWorkSpace[args];
936 AT.pWorkSpace[args] = arg1;
937 AT.pWorkSpace[args+1] = arg2;
939 m = mh = AT.WorkPointer;
940 mm = arg1+ARGHEAD; i = *mm;
961 for ( i = 0; i < numargs; i++ ) {
962 arg1 = AT.pWorkSpace[args+i];
964 m = (WORD *)Malloc1(*arg1*
sizeof(WORD),
"argbuffer type 0");
974 else if ( *arg1 == -DOLLAREXPRESSION ) {
975 d = DolToTerms(BHEAD arg1[1]);
976 abuf[i].buffer = d->where;
980 if ( *m ) argsdone++;
982 abuf[i].size = m-abuf[i].buffer;
984 else if ( *arg1 == -EXPRESSION ) {
985 abuf[i].buffer = CreateExpression(BHEAD arg1[1]);
988 if ( *m ) argsdone++;
990 abuf[i].size = m-abuf[i].buffer;
993 MLOCK(ErrorMessageLock);
994 MesPrint(
"What argument is this?");
995 MUNLOCK(ErrorMessageLock);
999 for ( i = 0; i < numargs; i++ ) {
1000 arg1 = abuf[i].buffer;
1001 if ( *arg1 == 0 ) {}
1002 else if ( arg1[*arg1] == 0 ) {
1006 ab = abuf[i]; abuf[i] = abuf[0]; abuf[0] = ab;
1007 mh = abuf[0].buffer;
1008 for ( j = 1; j < numargs; j++ ) {
1011 GCDterms(BHEAD mh,m,mh); m += *m;
1013 if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
1014 gcdisone = 1; sign = 1;
break;
1019 mm = mh + *mh;
if ( mm[-1] < 0 ) { sign = -1; mm[-1] = -mm[-1]; }
1020 mstop = mm - mm[-1]; m = mh+1; mlength = mm[-1];
1021 while ( tin < t ) *tout++ = *tin++;
1022 while ( m < mstop ) *tout++ = *m++;
1024 while ( tin < tstop ) *tout++ = *tin++;
1025 tlength = REDLENG(tlength);
1026 mlength = REDLENG(mlength);
1027 if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)mstop,mlength,
1028 (UWORD *)tout,&newlength) < 0 )
goto CalledFrom;
1029 mlength = INCLENG(newlength);
1030 tout += ABS(mlength);
1031 tout[-1] = mlength*sign;
1032 *termout = tout - termout;
1033 AT.WorkPointer = tout;
1034 if ( argsdone &&
Generator(BHEAD termout,level) < 0 )
goto CalledFrom;
1042 for ( i = 1; i < numargs; i++ ) {
1043 for ( ii = i; ii > 0; ii-- ) {
1044 if ( abuf[ii-1].size <= abuf[ii].size )
break;
1045 ab = abuf[ii-1]; abuf[ii-1] = abuf[ii]; abuf[ii] = ab;
1053 gcdout = abuf[ii].buffer;
1054 for ( i = 0; i < numargs; i++ ) {
1055 if ( abuf[i].buffer[0] ) { gcdout = abuf[i].buffer; ii = i; i++; argsdone++;
break; }
1057 for ( ; i < numargs; i++ ) {
1058 if ( abuf[i].buffer[0] ) {
1059 g = GCDfunction3(BHEAD gcdout,abuf[i].buffer);
1061 if ( gcdout != abuf[ii].buffer ) M_free(gcdout,
"gcdout");
1063 if ( gcdout[*gcdout] == 0 && gcdout[0] == 4 && gcdout[1] == 1
1064 && gcdout[2] == 1 && gcdout[3] == 3 )
break;
1069 tlength = REDLENG(tlength);
1071 tin = term; tout = termout;
while ( tin < t ) *tout++ = *tin++;
1073 mnext = mm + *mm; mlength = mnext[-1]; mstop = mnext - ABS(mlength);
1075 while ( mm < mstop ) *tout++ = *mm++;
1076 while ( tin < tstop ) *tout++ = *tin++;
1077 mlength = REDLENG(mlength);
1078 if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)mm,mlength,
1079 (UWORD *)tout,&newlength) < 0 )
goto CalledFrom;
1080 mlength = INCLENG(newlength);
1081 tout += ABS(mlength);
1083 *termout = tout - termout;
1084 AT.WorkPointer = tout;
1085 if ( argsdone &&
Generator(BHEAD termout,level) < 0 )
goto CalledFrom;
1088 if ( action && ( gcdout != abuf[ii].buffer ) ) M_free(gcdout,
"gcdout");
1095 for ( i = 0; i < numargs; i++ ) {
1096 if ( abuf[i].type == 0 ) { M_free(abuf[i].buffer,
"argbuffer type 0"); }
1097 else if ( abuf[i].type == 1 ) {
1099 if ( d->factors ) M_free(d->factors,
"Dollar factors");
1100 M_free(d,
"Copy of dollar variable");
1102 else if ( abuf[i].type == 2 ) { M_free(abuf[i].buffer,
"CreateExpression"); }
1104 M_free(abuf,
"argbuffer");
1109 AT.pWorkPointer = args;
1110 AT.WorkPointer = termout;
1114 MLOCK(ErrorMessageLock);
1115 MesCall(
"GCDfunction");
1116 MUNLOCK(ErrorMessageLock);
1141WORD *GCDfunction3(PHEAD WORD *in1, WORD *in2)
1144 WORD oldsorttype = AR.SortType, *ow = AT.WorkPointer;;
1145 WORD *t, *tt, *gcdout, *term1, *term2, *confree1, *confree2, *gcdout1, *proper1, *proper2;
1146 int i, actionflag1, actionflag2;
1147 WORD startebuf = cbuf[AT.ebufnum].numrhs;
1148 WORD tryterm1, tryterm2;
1149 if ( in2[*in2] == 0 ) { t = in1; in1 = in2; in2 = t; }
1150 if ( in1[*in1] == 0 ) {
1151 gcdout = (WORD *)Malloc1((*in1+1)*
sizeof(WORD),
"gcdout");
1152 i = *in1; t = gcdout; tt = in1; NCOPY(t,tt,i); *t = 0;
1155 GCDterms(BHEAD gcdout,t,gcdout);
1156 if ( gcdout[0] == 4 && gcdout[1] == 1
1157 && gcdout[2] == 1 && gcdout[3] == 3 )
break;
1160 AT.WorkPointer = ow;
1167 AR.SortType = SORTHIGHFIRST;
1168 term1 = TermMalloc(
"GCDfunction3-a");
1169 term2 = TermMalloc(
"GCDfunction3-b");
1172 tryterm1 = AN.tryterm; AN.tryterm = 0;
1174 tryterm2 = AN.tryterm; AN.tryterm = 0;
1179 GCDterms(BHEAD term1,term2,term1);
1180 TermFree(term2,
"GCDfunction3-b");
1185 if ( ( proper1 = PutExtraSymbols(BHEAD confree1,startebuf,&actionflag1) ) == 0 )
goto CalledFrom;
1186 if ( confree1 != in1 ) {
1187 if ( tryterm1 ) { TermFree(confree1,
"TakeContent"); }
1188 else { M_free(confree1,
"TakeContent"); }
1193 if ( ( proper2 = PutExtraSymbols(BHEAD confree2,startebuf,&actionflag2) ) == 0 )
goto CalledFrom;
1194 if ( confree2 != in2 ) {
1195 if ( tryterm2 ) { TermFree(confree2,
"TakeContent"); }
1196 else { M_free(confree2,
"TakeContent"); }
1204 gcdout1 =
poly_gcd(BHEAD proper1,proper2,0);
1205 M_free(proper1,
"PutExtraSymbols");
1206 M_free(proper2,
"PutExtraSymbols");
1208 AR.SortType = oldsorttype;
1209 if ( actionflag1 || actionflag2 ) {
1210 if ( ( gcdout = TakeExtraSymbols(BHEAD gcdout1,startebuf) ) == 0 )
goto CalledFrom;
1211 M_free(gcdout1,
"gcdout");
1217 cbuf[AT.ebufnum].numrhs = startebuf;
1221 if ( term1[0] != 4 || term1[3] != 3 || term1[1] != 1 || term1[2] != 1 ) {
1223 if ( ( gcdout1 = MultiplyWithTerm(BHEAD gcdout,term1,2) ) == 0 )
goto CalledFrom;
1225 M_free(gcdout,
"gcdout");
1228 TermFree(term1,
"GCDfunction3-a");
1229 AT.WorkPointer = ow;
1233 MLOCK(ErrorMessageLock);
1234 MesCall(
"GCDfunction3");
1235 MUNLOCK(ErrorMessageLock);
1244WORD *PutExtraSymbols(PHEAD WORD *in,WORD startebuf,
int *actionflag)
1246 WORD *termout = AT.WorkPointer;
1255 if ( action > 0 ) *actionflag = 1;
1259 if (
EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 )
goto CalledFrom;
1262 MLOCK(ErrorMessageLock);
1263 MesCall(
"PutExtraSymbols");
1264 MUNLOCK(ErrorMessageLock);
1273WORD *TakeExtraSymbols(PHEAD WORD *in,WORD startebuf)
1275 CBUF *C = cbuf+AC.cbufnum;
1276 CBUF *CC = cbuf+AT.ebufnum;
1277 WORD *oldworkpointer = AT.WorkPointer, *termout;
1279 termout = AT.WorkPointer;
1282 if ( ConvertFromPoly(BHEAD in,termout,numxsymbol,CC->numrhs-startebuf+numxsymbol,startebuf-numxsymbol,1) <= 0 ) {
1287 AT.WorkPointer = termout + *termout;
1291 if (
Generator(BHEAD termout,C->numlhs) ) {
1296 AT.WorkPointer = oldworkpointer;
1297 if (
EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 )
goto CalledFrom;
1301 MLOCK(ErrorMessageLock);
1302 MesCall(
"TakeExtraSymbols");
1303 MUNLOCK(ErrorMessageLock);
1312WORD *MultiplyWithTerm(PHEAD WORD *in, WORD *term, WORD par)
1314 WORD *termout, *t, *tt, *tstop, *ttstop;
1315 WORD length, length1, length2;
1316 WORD oldsorttype = AR.SortType;
1317 COMPARE oldcompareroutine = AR.CompareRoutine;
1320 if ( par == 0 || par == 2 ) AR.SortType = SORTHIGHFIRST;
1321 else AR.SortType = SORTLOWFIRST;
1322 termout = AT.WorkPointer;
1326 tstop = in + *in; tstop -= ABS(tstop[-1]); t = in + 1;
1327 while ( t < tstop ) *tt++ = *t++;
1328 ttstop = term + *term; ttstop -= ABS(ttstop[-1]); t = term + 1;
1329 while ( t < ttstop ) *tt++ = *t++;
1330 length1 = REDLENG(in[*in-1]); length2 = REDLENG(term[*term-1]);
1331 if ( MulRat(BHEAD (UWORD *)tstop,length1,
1332 (UWORD *)ttstop,length2,(UWORD *)tt,&length) )
goto CalledFrom;
1333 length = INCLENG(length);
1334 tt += ABS(length); tt[-1] = length;
1335 *termout = tt - termout;
1343 if (
EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 )
goto CalledFrom;
1346 if (
EndSort(BHEAD termout,1) < 0 )
goto CalledFrom;
1349 AR.CompareRoutine = oldcompareroutine;
1351 AR.SortType = oldsorttype;
1355 MLOCK(ErrorMessageLock);
1356 MesCall(
"MultiplyWithTerm");
1357 MUNLOCK(ErrorMessageLock);
1379 WORD *t, *tstop, *tcom, *tout, *tstore, *r, *rstop, *m, *mm, *w, *ww, *wterm;
1380 WORD *tnext, *tt, *tterm, code[2];
1382 int i, j, k, action = 0, sign;
1383 UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMbuffer2, *ap;
1384 WORD GCDlen, GCDlen2, LCMlen, LCMlen2, length, redlength, len1, len2;
1385 tout = tstore = term+1;
1391 tstop = tnext-ABS(tnext[-1]);
1393 while ( t < tstop ) {
1394 if ( *t == INDEX ) {
1395 i = t[1]; NCOPY(tout,t,i);
break;
1399 if ( tout > tstore ) {
1403 rstop = tnext - ABS(tnext[-1]);
1405 if ( r == rstop )
goto noindices;
1406 while ( r < rstop ) {
1407 if ( *r != INDEX ) { r += r[1];
continue; }
1409 while ( m < tout ) {
1410 for ( i = 2; i < r[1]; i++ ) {
1411 if ( *m == r[i] )
break;
1412 if ( *m > r[i] )
continue;
1414 while ( mm < tout ) { mm[-1] = mm[0]; mm++; }
1415 tout--; tstore[1]--; m--;
1421 if ( r >= rstop || tout <= tstore+2 ) {
1422 tout = tstore;
break;
1425 if ( tout > tstore+2 ) {
1429 tnext = t + *t; t++; w++;
1430 while ( *t != INDEX ) { i = t[1]; NCOPY(w,t,i); }
1431 tt = t + t[1]; t += 2; r = tstore+2; ww = w; *w++ = INDEX; w++;
1432 while ( r < tout && t < tt ) {
1433 if ( *r > *t ) { *w++ = *t++; }
1434 else if ( *r == *t ) { r++; t++; }
1435 else goto CalledFrom;
1437 if ( r < tout )
goto CalledFrom;
1438 while ( t < tt ) *w++ = *t++;
1440 if ( ww[1] == 2 ) w = ww;
1441 while ( t < tnext ) *w++ = *t++;
1453 code[0] = VECTOR; code[1] = DELTA;
1454 for ( k = 0; k < 2; k++ ) {
1457 tstop = tnext-ABS(tnext[-1]);
1459 while ( t < tstop ) {
1460 if ( *t == code[k] ) {
1461 i = t[1]; NCOPY(tout,t,i);
break;
1465 if ( tout > tstore ) {
1469 rstop = tnext - ABS(tnext[-1]);
1471 if ( r == rstop ) { tstore = tout;
goto novectors; }
1472 while ( r < rstop ) {
1473 if ( *r != code[k] ) { r += r[1];
continue; }
1475 while ( m < tout ) {
1476 for ( i = 2; i < r[1]; i += 2 ) {
1477 if ( *m == r[i] && m[1] == r[i+1] )
break;
1478 if ( *m > r[i] || ( *m == r[i] && m[1] > r[i+1] ) )
continue;
1480 while ( mm < tout ) { mm[-2] = mm[0]; mm[-1] = mm[1]; mm += 2; }
1481 tout -= 2; tstore[1] -= 2; m -= 2;
1487 if ( r >= rstop || tout <= tstore+2 ) {
1488 tout = tstore;
break;
1491 if ( tout > tstore+2 ) {
1495 tnext = t + *t; t++; w++;
1496 while ( *t != code[k] ) { i = t[1]; NCOPY(w,t,i); }
1497 tt = t + t[1]; t += 2; r = tstore+2; ww = w; *w++ = code[k]; w++;
1498 while ( r < tout && t < tt ) {
1499 if ( ( *r > *t ) || ( *r == *t && r[1] > t[1] ) )
1500 { *w++ = *t++; *w++ = *t++; }
1501 else if ( *r == *t && r[1] == t[1] ) { r += 2; t += 2; }
1502 else goto CalledFrom;
1504 if ( r < tout )
goto CalledFrom;
1505 while ( t < tt ) *w++ = *t++;
1507 if ( ww[1] == 2 ) w = ww;
1508 while ( t < tnext ) *w++ = *t++;
1523 tstop = tnext-ABS(tnext[-1]);
1526 while ( t < tstop ) {
1527 if ( *t >= FUNCTION ) {
1528 if ( functions[*t-FUNCTION].commute ) {
1529 if ( tcom == 0 ) { tcom = tstore; }
1531 for ( i = 0; i < t[1]; i++ ) {
1532 if ( t[i] != tcom[i] ) {
1533 MLOCK(ErrorMessageLock);
1534 MesPrint(
"GCD or factorization of more than one noncommuting object not allowed");
1535 MUNLOCK(ErrorMessageLock);
1541 i = t[1]; NCOPY(tstore,t,i);
1545 if ( tout > tstore ) {
1548 tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1549 if ( t == tstop )
goto nofunctions;
1551 while ( r < tout ) {
1553 while ( tt < tstop ) {
1554 for ( i = 0; i < r[1]; i++ ) {
1555 if ( r[i] != tt[i] )
break;
1557 if ( i == r[1] ) { r += r[1];
goto nextr1; }
1562 m = r; mm = r + r[1];
1563 while ( mm < tout ) *m++ = *mm++;
1567 if ( tout <= tstore )
break;
1571 if ( tout > tstore ) {
1577 while ( r < tout ) {
1578 t = in; ww = in; w = ww+1;
1583 for ( i = 0; i < r[1]; i++ ) {
1584 if ( t[i] != r[i] ) {
1585 j = t[1]; NCOPY(w,t,j);
1591 while ( t < tnext ) *w++ = *t++;
1612 tterm = AT.WorkPointer; tt = tterm+1;
1613 tout[0] = SYMBOL; tout[1] = 2;
1615 tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1616 while ( t < tstop ) {
1617 if ( *t == SYMBOL ) {
1618 for ( i = 0; i < t[1]; i++ ) tout[i] = t[i];
1625 tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1631 while ( t < tstop ) {
1632 if ( *t == SYMBOL ) {
1633 MergeSymbolLists(BHEAD tout,t,-1);
1641 if ( tout[1] > 2 ) {
1643 tt[0] = t[0]; tt[1] = t[1];
1644 for ( i = 2; i < t[1]; i += 2 ) {
1645 tt[i] = t[i]; tt[i+1] = -t[i+1];
1660 tout[0] = DOTPRODUCT; tout[1] = 2;
1663 tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1668 while ( t < tstop ) {
1669 if ( *t == DOTPRODUCT ) {
1670 MergeDotproductLists(BHEAD tout,t,-1);
1677 if ( tout[1] > 2 ) {
1679 tt[0] = t[0]; tt[1] = t[1];
1680 for ( i = 2; i < t[1]; i += 3 ) {
1681 tt[i] = t[i]; tt[i+1] = t[i+1]; tt[i+2] = -t[i+2];
1694 AT.WorkPointer = tt;
1695 if ( AN.cmod != 0 ) {
1697 t = in; tnext = t + *t; tstop = tnext - ABS(tnext[-1]);
1699 if ( tnext[-1] < 0 ) x += AC.cmod[0];
1700 if (
GetModInverses(x,(WORD)(AN.cmod[0]),&ix,&ip) )
goto CalledFrom;
1701 *tout++ = x; *tout++ = 1; *tout++ = 3;
1702 *tt++ = ix; *tt++ = 1; *tt++ = 3;
1705 GCDbuffer = NumberMalloc(
"MakeInteger");
1706 GCDbuffer2 = NumberMalloc(
"MakeInteger");
1707 LCMbuffer = NumberMalloc(
"MakeInteger");
1708 LCMbuffer2 = NumberMalloc(
"MakeInteger");
1710 tnext = t + *t; length = tnext[-1];
1711 if ( length < 0 ) { sign = -1; length = -length; }
1713 tstop = tnext - length;
1714 redlength = (length-1)/2;
1715 for ( i = 0; i < redlength; i++ ) {
1716 GCDbuffer[i] = (UWORD)(tstop[i]);
1717 LCMbuffer[i] = (UWORD)(tstop[redlength+i]);
1719 GCDlen = LCMlen = redlength;
1720 while ( GCDbuffer[GCDlen-1] == 0 ) GCDlen--;
1721 while ( LCMbuffer[LCMlen-1] == 0 ) LCMlen--;
1724 tnext = t + *t; length = ABS(tnext[-1]);
1725 tstop = tnext - length; redlength = (length-1)/2;
1726 len1 = len2 = redlength;
1727 den = tstop + redlength;
1728 while ( tstop[len1-1] == 0 ) len1--;
1729 while ( den[len2-1] == 0 ) len2--;
1730 if ( GCDlen == 1 && GCDbuffer[0] == 1 ) {}
1732 GcdLong(BHEAD (UWORD *)tstop,len1,GCDbuffer,GCDlen,GCDbuffer2,&GCDlen2);
1733 ap = GCDbuffer; GCDbuffer = GCDbuffer2; GCDbuffer2 = ap;
1734 a = GCDlen; GCDlen = GCDlen2; GCDlen2 = a;
1736 if ( len2 == 1 && den[0] == 1 ) {}
1738 GcdLong(BHEAD LCMbuffer,LCMlen,(UWORD *)den,len2,LCMbuffer2,&LCMlen2);
1739 DivLong((UWORD *)den,len2,LCMbuffer2,LCMlen2,
1740 GCDbuffer2,&GCDlen2,(UWORD *)AT.WorkPointer,&a);
1741 MulLong(LCMbuffer,LCMlen,GCDbuffer2,GCDlen2,LCMbuffer2,&LCMlen2);
1742 ap = LCMbuffer; LCMbuffer = LCMbuffer2; LCMbuffer2 = ap;
1743 a = LCMlen; LCMlen = LCMlen2; LCMlen2 = a;
1747 if ( GCDlen != 1 || GCDbuffer[0] != 1 || LCMlen != 1 || LCMbuffer[0] != 1 ) {
1748 redlength = GCDlen;
if ( LCMlen > GCDlen ) redlength = LCMlen;
1749 for ( i = 0; i < GCDlen; i++ ) *tout++ = (WORD)(GCDbuffer[i]);
1750 for ( ; i < redlength; i++ ) *tout++ = 0;
1751 for ( i = 0; i < LCMlen; i++ ) *tout++ = (WORD)(LCMbuffer[i]);
1752 for ( ; i < redlength; i++ ) *tout++ = 0;
1753 *tout++ = (2*redlength+1)*sign;
1754 for ( i = 0; i < LCMlen; i++ ) *tt++ = (WORD)(LCMbuffer[i]);
1755 for ( ; i < redlength; i++ ) *tt++ = 0;
1756 for ( i = 0; i < GCDlen; i++ ) *tt++ = (WORD)(GCDbuffer[i]);
1757 for ( ; i < redlength; i++ ) *tt++ = 0;
1758 *tt++ = (2*redlength+1)*sign;
1762 *tout++ = 1; *tout++ = 1; *tout++ = 3*sign;
1763 *tt++ = 1; *tt++ = 1; *tt++ = 3*sign;
1764 if ( sign != 1 ) action++;
1767 NumberFree(LCMbuffer2,
"MakeInteger");
1768 NumberFree(LCMbuffer ,
"MakeInteger");
1769 NumberFree(GCDbuffer2,
"MakeInteger");
1770 NumberFree(GCDbuffer ,
"MakeInteger");
1777 *tterm = tt - tterm;
1778 AT.WorkPointer = tt;
1779 inp = MultiplyWithTerm(BHEAD in,tterm,2);
1780 AT.WorkPointer = tterm;
1786 *term = tout - term;
1787 AT.WorkPointer = tterm;
1790 MLOCK(ErrorMessageLock);
1791 MesCall(
"TakeContent");
1792 MUNLOCK(ErrorMessageLock);
1808int MergeSymbolLists(PHEAD WORD *old, WORD *extra,
int par)
1811 WORD *
new = TermMalloc(
"MergeSymbolLists");
1812 WORD *t1, *t2, *fill;
1814 fill =
new + 2; *
new = SYMBOL;
1815 i1 = old[1] - 2; i2 = extra[1] - 2;
1816 t1 = old + 2; t2 = extra + 2;
1819 while ( i1 > 0 && i2 > 0 ) {
1821 if ( t2[1] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1825 else if ( *t1 < *t2 ) {
1826 if ( t1[1] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1830 else if ( t1[1] < t2[1] ) {
1831 *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1835 *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1839 for ( ; i1 > 0; i1 -= 2 ) {
1840 if ( t1[1] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1843 for ( ; i2 > 0; i2 -= 2 ) {
1844 if ( t2[1] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1849 while ( i1 > 0 && i2 > 0 ) {
1851 if ( t2[1] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1855 else if ( *t1 < *t2 ) {
1856 if ( t1[1] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1860 else if ( t1[1] > t2[1] ) {
1861 *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1865 *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1869 for ( ; i1 > 0; i1 -= 2 ) {
1870 if ( t1[1] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1873 for ( ; i2 > 0; i2 -= 2 ) {
1874 if ( t2[1] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1879 while ( i1 > 0 && i2 > 0 ) {
1883 else if ( *t1 < *t2 ) {
1886 else if ( ( t1[1] > 0 ) && ( t2[1] < 0 ) ) { t1 += 2; t2 += 2; i1 -= 2; i2 -= 2; }
1887 else if ( ( t1[1] < 0 ) && ( t2[1] > 0 ) ) { t1 += 2; t2 += 2; i1 -= 2; i2 -= 2; }
1888 else if ( t1[1] > 0 ) {
1889 if ( t1[1] < t2[1] ) {
1890 *fill++ = *t1++; *fill++ = *t1++; t2 += 2; i2 -= 2;
1893 *fill++ = *t2++; *fill++ = *t2++; t1 += 2; i1 -= 2;
1897 if ( t2[1] < t1[1] ) {
1898 *fill++ = *t2++; *fill++ = *t2++; t1 += 2; i1 -= 2; i2 -= 2;
1901 *fill++ = *t1++; *fill++ = *t1++; t2 += 2; i1 -= 2; i2 -= 2;
1905 for ( ; i1 > 0; i1-- ) *fill++ = *t1++;
1906 for ( ; i2 > 0; i2-- ) *fill++ = *t2++;
1909 i1 =
new[1] = fill -
new;
1910 t2 =
new; t1 = old; NCOPY(t1,t2,i1);
1911 TermFree(
new,
"MergeSymbolLists");
1927int MergeDotproductLists(PHEAD WORD *old, WORD *extra,
int par)
1930 WORD *
new = TermMalloc(
"MergeDotproductLists");
1931 WORD *t1, *t2, *fill;
1934 i1 = old[1] - 2; i2 = extra[1] - 2;
1935 t1 = old + 2; t2 = extra + 2;
1938 while ( i1 > 0 && i2 > 0 ) {
1939 if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1940 if ( t2[2] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; }
1943 else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1944 if ( t1[2] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; }
1947 else if ( t1[2] < t2[2] ) {
1948 *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1951 *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1956 while ( i1 > 0 && i2 > 0 ) {
1957 if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1958 if ( t2[2] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; }
1961 else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1962 if ( t1[2] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; }
1965 else if ( t1[2] > t2[2] ) {
1966 *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1969 *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1974 while ( i1 > 0 && i2 > 0 ) {
1975 if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1978 else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1981 else if ( ( t1[2] > 0 ) && ( t2[2] < 0 ) ) { t1 += 3; t2 += 3; }
1982 else if ( ( t1[2] < 0 ) && ( t2[2] > 0 ) ) { t1 += 3; t2 += 3; }
1983 else if ( t1[2] > 0 ) {
1984 if ( t1[2] < t2[2] ) {
1985 *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1988 *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1992 if ( t2[2] < t1[2] ) {
1993 *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1996 *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
2002 i1 =
new[1] = fill -
new;
2003 t2 =
new; t1 = old; NCOPY(t1,t2,i1);
2004 TermFree(
new,
"MergeDotproductLists");
2020WORD *CreateExpression(PHEAD WORD nexp)
2023 CBUF *C = cbuf+AC.cbufnum;
2024 POSITION startposition, oldposition;
2026 WORD *term, *oldipointer = AR.CompressPointer;
2028 switch ( Expressions[nexp].status ) {
2029 case HIDDENLEXPRESSION:
2030 case HIDDENGEXPRESSION:
2031 case DROPHLEXPRESSION:
2032 case DROPHGEXPRESSION:
2033 case UNHIDELEXPRESSION:
2034 case UNHIDEGEXPRESSION:
2035 AR.GetOneFile = 2; fi = AR.hidefile;
2038 AR.GetOneFile = 0; fi = AR.infile;
2041 SeekScratch(fi,&oldposition);
2042 startposition = AS.OldOnFile[nexp];
2043 term = AT.WorkPointer;
2044 if ( GetOneTerm(BHEAD term,fi,&startposition,0) <= 0 )
goto CalledFrom;
2046 AR.CompressPointer = oldipointer;
2047 while ( GetOneTerm(BHEAD term,fi,&startposition,0) > 0 ) {
2048 AT.WorkPointer = term + *term;
2049 if (
Generator(BHEAD term,C->numlhs) ) {
2053 AR.CompressPointer = oldipointer;
2055 AT.WorkPointer = term;
2056 if (
EndSort(BHEAD (WORD *)((VOID *)(&term)),2) < 0 )
goto CalledFrom;
2057 SetScratch(fi,&oldposition);
2060 MLOCK(ErrorMessageLock);
2061 MesCall(
"CreateExpression");
2062 MUNLOCK(ErrorMessageLock);
2076int GCDterms(PHEAD WORD *term1, WORD *term2, WORD *termout)
2079 WORD *t1, *t1stop, *t1next, *t2, *t2stop, *t2next, *tout, *tt1, *tt2;
2080 int count1, count2, i, ii, x1, sign;
2081 WORD length1, length2;
2082 t1 = term1 + *term1; t1stop = t1 - ABS(t1[-1]); t1 = term1+1;
2083 t2 = term2 + *term2; t2stop = t2 - ABS(t2[-1]); t2 = term2+1;
2085 while ( t1 < t1stop ) {
2086 t1next = t1 + t1[1];
2088 if ( *t1 == SYMBOL ) {
2089 while ( t2 < t2stop && *t2 != SYMBOL ) t2 += t2[1];
2090 if ( t2 < t2stop && *t2 == SYMBOL ) {
2092 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2093 while ( tt1 < t1next && tt2 < t2next ) {
2094 if ( *tt1 < *tt2 ) tt1 += 2;
2095 else if ( *tt1 > *tt2 ) tt2 += 2;
2096 else if ( ( tt1[1] > 0 && tt2[1] < 0 ) ||
2097 ( tt2[1] > 0 && tt1[1] < 0 ) ) {
2102 if ( tt1[1] < 0 ) {
if ( tt2[1] > x1 ) x1 = tt2[1]; }
2103 else {
if ( tt2[1] < x1 ) x1 = tt2[1]; }
2104 tout[count1+2] = *tt1;
2105 tout[count1+3] = x1;
2111 *tout = SYMBOL; tout[1] = count1+2; tout += tout[1];
2115 else if ( *t1 == DOTPRODUCT ) {
2116 while ( t2 < t2stop && *t2 != DOTPRODUCT ) t2 += t2[1];
2117 if ( t2 < t2stop && *t2 == DOTPRODUCT ) {
2119 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2120 while ( tt1 < t1next && tt2 < t2next ) {
2121 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 3;
2122 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 3;
2123 else if ( ( tt1[2] > 0 && tt2[2] < 0 ) ||
2124 ( tt2[2] > 0 && tt1[2] < 0 ) ) {
2129 if ( tt1[2] < 0 ) {
if ( tt2[2] > x1 ) x1 = tt2[2]; }
2130 else {
if ( tt2[2] < x1 ) x1 = tt2[2]; }
2131 tout[count1+2] = *tt1;
2132 tout[count1+3] = tt1[1];
2133 tout[count1+4] = x1;
2139 *tout = DOTPRODUCT; tout[1] = count1+2; tout += tout[1];
2143 else if ( *t1 == VECTOR ) {
2144 while ( t2 < t2stop && *t2 != VECTOR ) t2 += t2[1];
2145 if ( t2 < t2stop && *t2 == VECTOR ) {
2147 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2148 while ( tt1 < t1next && tt2 < t2next ) {
2149 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2150 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2152 tout[count1+2] = *tt1;
2153 tout[count1+3] = tt1[1];
2159 *tout = VECTOR; tout[1] = count1+2; tout += tout[1];
2163 else if ( *t1 == INDEX ) {
2164 while ( t2 < t2stop && *t2 != INDEX ) t2 += t2[1];
2165 if ( t2 < t2stop && *t2 == INDEX ) {
2167 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2168 while ( tt1 < t1next && tt2 < t2next ) {
2169 if ( *tt1 < *tt2 ) tt1 += 1;
2170 else if ( *tt1 > *tt2 ) tt2 += 1;
2172 tout[count1+2] = *tt1;
2178 *tout = INDEX; tout[1] = count1+2; tout += tout[1];
2182 else if ( *t1 == DELTA ) {
2183 while ( t2 < t2stop && *t2 != DELTA ) t2 += t2[1];
2184 if ( t2 < t2stop && *t2 == DELTA ) {
2186 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2187 while ( tt1 < t1next && tt2 < t2next ) {
2188 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2189 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2191 tout[count1+2] = *tt1;
2192 tout[count1+3] = tt1[1];
2198 *tout = DELTA; tout[1] = count1+2; tout += tout[1];
2202 else if ( *t1 >= FUNCTION ) {
2208 while ( t1next < t1stop && *t1 == *t1next && t1[1] == t1next[1] ) {
2209 for ( i = 2; i < t1[1]; i++ ) {
2210 if ( t1[i] != t1next[i] )
break;
2212 if ( i < t1[1] )
break;
2214 t1next += t1next[1];
2217 while ( t2 < t2stop ) {
2218 if ( *t2 == *t1 && t2[1] == t1[1] ) {
2219 for ( i = 2; i < t1[1]; i++ ) {
2220 if ( t2[i] != t1[i] )
break;
2222 if ( i >= t1[1] ) count2++;
2226 if ( count1 < count2 ) count2 = count1;
2229 while ( count2 > 0 ) { tout += tout[1]; count2--; }
2243 length1 = term1[*term1-1]; ii = i = ABS(length1); t1 = t1stop;
2244 if ( t1 != tout ) { NCOPY(tout,t1,i); tout -= ii; }
2245 length2 = term2[*term2-1];
2246 if ( length1 < 0 && length2 < 0 ) sign = -1;
2247 if ( AccumGCD(BHEAD (UWORD *)tout,&length1,(UWORD *)t2stop,length2) ) {
2248 MLOCK(ErrorMessageLock);
2249 MesCall(
"GCDterms");
2250 MUNLOCK(ErrorMessageLock);
2253 if ( sign < 0 && length1 > 0 ) length1 = -length1;
2254 tout += ABS(length1); tout[-1] = length1;
2255 *termout = tout - termout; *tout = 0;
2264int ReadPolyRatFun(PHEAD WORD *term)
2266 WORD *oldworkpointer = AT.WorkPointer;
2268 WORD *t, *fun, *nextt, *num, *den, *t1, *t2, size, numsize, densize;
2269 WORD *term1, *term2, *confree1, *confree2, *gcd, *num1, *den1, move, *newnum, *newden;
2270 WORD *tstop, *m1, *m2;
2271 WORD oldsorttype = AR.SortType;
2272 COMPARE oldcompareroutine = AR.CompareRoutine;
2273 AR.SortType = SORTHIGHFIRST;
2276 tstop = term + *term; tstop -= ABS(tstop[-1]);
2277 if ( term + *term == AT.WorkPointer ) flag = 1;
2280 while ( t < tstop ) {
2281 if ( *t != AR.PolyFun ) { t += t[1];
continue; }
2282 if ( ( t[2] & MUSTCLEANPRF ) == 0 ) { t += t[1];
continue; }
2285 if ( fun[1] > FUNHEAD && fun[FUNHEAD] == -SNUMBER && fun[FUNHEAD+1] == 0 )
2286 { *term = 0;
break; }
2287 if ( FromPolyRatFun(BHEAD fun, &num, &den) > 0 ) { t = nextt;
continue; }
2288 if ( *num == ARGHEAD ) { *term = 0;
break; }
2295 term1 = TermMalloc(
"ReadPolyRatFun");
2296 term2 = TermMalloc(
"ReadPolyRatFun");
2299 GCDclean(BHEAD term1,term2);
2301 gcd =
poly_gcd(BHEAD confree1,confree2,1);
2302 newnum = PolyDiv(BHEAD confree1,gcd,
"ReadPolyRatFun");
2303 newden = PolyDiv(BHEAD confree2,gcd,
"ReadPolyRatFun");
2304 TermFree(confree2,
"ReadPolyRatFun");
2305 TermFree(confree1,
"ReadPolyRatFun");
2306 num1 = MULfunc(BHEAD term1,newnum);
2307 den1 = MULfunc(BHEAD term2,newden);
2308 TermFree(newnum,
"ReadPolyRatFun");
2309 TermFree(newden,
"ReadPolyRatFun");
2311 TermFree(gcd,
"poly_gcd");
2312 TermFree(term1,
"ReadPolyRatFun");
2313 TermFree(term2,
"ReadPolyRatFun");
2320 if ( num1[0] == 4 && num1[4] == 0 && num1[2] == 1 && num1[1] > 0 ) {
2321 numsize = 2; num1[0] = -SNUMBER;
2322 if ( num1[3] < 0 ) num1[1] = -num1[1];
2324 else if ( num1[0] == 8 && num1[8] == 0 && num1[7] == 3 && num1[6] == 1
2325 && num1[5] == 1 && num1[1] == SYMBOL && num1[4] == 1 ) {
2326 numsize = 2; num1[0] = -SYMBOL; num1[1] = num1[3];
2328 else { m1 = num1;
while ( *m1 ) m1 += *m1; numsize = (m1-num1)+ARGHEAD; }
2329 if ( den1[0] == 4 && den1[4] == 0 && den1[2] == 1 && den1[1] > 0 ) {
2330 densize = 2; den1[0] = -SNUMBER;
2331 if ( den1[3] < 0 ) den1[1] = -den1[1];
2333 else if ( den1[0] == 8 && den1[8] == 0 && den1[7] == 3 && den1[6] == 1
2334 && den1[5] == 1 && den1[1] == SYMBOL && den1[4] == 1 ) {
2335 densize = 2; den1[0] = -SYMBOL; den1[1] = den1[3];
2337 else { m2 = den1;
while ( *m2 ) m2 += *m2; densize = (m2-den1)+ARGHEAD; }
2338 size = FUNHEAD+numsize+densize;
2340 if ( size > fun[1] ) {
2341 move = size - fun[1];
2342 t1 = term+*term; t2 = t1+move;
2343 while ( t1 > nextt ) *--t2 = *--t1;
2344 tstop += move; nextt += move;
2347 else if ( size < fun[1] ) {
2349 t2 = fun+size; t1 = nextt;
2350 tstop -= move; nextt -= move;
2352 while ( t1 < t ) *t2++ = *t1++;
2356 fun[1] = size; fun[2] = 0;
2357 t2 = fun+FUNHEAD; t1 = num1;
2358 if ( *num1 < 0 ) { *t2++ = num1[0]; *t2++ = num1[1]; }
2359 else { *t2++ = numsize; *t2++ = 0; FILLARG(t2);
2360 i = numsize-ARGHEAD; NCOPY(t2,t1,i) }
2362 if ( *den1 < 0 ) { *t2++ = den1[0]; *t2++ = den1[1]; }
2363 else { *t2++ = densize; *t2++ = 0; FILLARG(t2);
2364 i = densize-ARGHEAD; NCOPY(t2,t1,i) }
2366 TermFree(num1,
"MULfunc");
2367 TermFree(den1,
"MULfunc");
2371 if ( flag ) AT.WorkPointer = term +*term;
2372 else AT.WorkPointer = oldworkpointer;
2373 AR.CompareRoutine = oldcompareroutine;
2374 AR.SortType = oldsorttype;
2383int FromPolyRatFun(PHEAD WORD *fun, WORD **numout, WORD **denout)
2385 WORD *nextfun, *tt, *num, *den;
2387 nextfun = fun + fun[1];
2389 num = AT.WorkPointer;
2391 if ( *fun != -SNUMBER && *fun != -SYMBOL )
goto Improper;
2392 ToGeneral(fun,num,0);
2393 tt = num + *num; *tt++ = 0;
2396 else { i = *fun; tt = num; NCOPY(tt,fun,i); *tt++ = 0; }
2399 if ( *fun != -SNUMBER && *fun != -SYMBOL )
goto Improper;
2400 ToGeneral(fun,den,0);
2401 tt = den + *den; *tt++ = 0;
2404 else { i = *fun; tt = den; NCOPY(tt,fun,i); *tt++ = 0; }
2405 *numout = num; *denout = den;
2406 if ( fun != nextfun ) {
return(1); }
2407 AT.WorkPointer = tt;
2410 MLOCK(ErrorMessageLock);
2411 MesPrint(
"Improper use of PolyRatFun");
2412 MesCall(
"FromPolyRatFun");
2413 MUNLOCK(ErrorMessageLock);
2437 WORD *t, *tstop, *tout, *tstore;
2438 WORD *tnext, *tt, *tterm;
2439 WORD *inp, a, *den, *oldworkpointer = AT.WorkPointer;
2440 int i, action = 0, sign, first;
2441 UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMbuffer2, *ap;
2442 WORD GCDlen, GCDlen2, LCMlen, LCMlen2, length, redlength, len1, len2;
2444 tout = tstore = term+1;
2453 tterm = AT.WorkPointer; tt = tterm+1;
2454 tout[0] = SYMBOL; tout[1] = 2;
2457 tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
2458 while ( t < tstop ) {
2460 if ( *t == SYMBOL ) {
2461 for ( i = 0; i < t[1]; i++ ) tout[i] = t[i];
2465 MLOCK(ErrorMessageLock);
2466 MesPrint ((
char*)
"ERROR: polynomials and polyratfuns must contain symbols only");
2467 MUNLOCK(ErrorMessageLock);
2471 else if ( *t == SYMBOL ) {
2472 MergeSymbolLists(BHEAD tout,t,-1);
2484 for ( i = 2; i < tout[1]; i += 2 ) {
2485 if ( tout[i+1] < 0 ) {
2486 if ( i == j ) { j += 2; }
2487 else { tout[j] = tout[i]; tout[j+1] = tout[i+1]; j += 2; }
2496 if ( tout[1] > 2 ) {
2498 tt[0] = t[0]; tt[1] = t[1];
2499 for ( i = 2; i < t[1]; i += 2 ) {
2500 tt[i] = t[i]; tt[i+1] = -t[i+1];
2513 AT.WorkPointer = tt;
2514 if ( AN.cmod != 0 ) {
2516 t = in; tnext = t + *t; tstop = tnext - ABS(tnext[-1]);
2518 if ( tnext[-1] < 0 ) x += AC.cmod[0];
2519 if (
GetModInverses(x,(WORD)(AN.cmod[0]),&ix,&ip) )
goto CalledFrom;
2520 *tout++ = x; *tout++ = 1; *tout++ = 3;
2521 *tt++ = ix; *tt++ = 1; *tt++ = 3;
2524 GCDbuffer = NumberMalloc(
"MakeInteger");
2525 GCDbuffer2 = NumberMalloc(
"MakeInteger");
2526 LCMbuffer = NumberMalloc(
"MakeInteger");
2527 LCMbuffer2 = NumberMalloc(
"MakeInteger");
2529 tnext = t + *t; length = tnext[-1];
2530 if ( length < 0 ) { sign = -1; length = -length; }
2532 tstop = tnext - length;
2533 redlength = (length-1)/2;
2534 for ( i = 0; i < redlength; i++ ) {
2535 GCDbuffer[i] = (UWORD)(tstop[i]);
2536 LCMbuffer[i] = (UWORD)(tstop[redlength+i]);
2538 GCDlen = LCMlen = redlength;
2539 while ( GCDbuffer[GCDlen-1] == 0 ) GCDlen--;
2540 while ( LCMbuffer[LCMlen-1] == 0 ) LCMlen--;
2543 tnext = t + *t; length = ABS(tnext[-1]);
2544 tstop = tnext - length; redlength = (length-1)/2;
2545 len1 = len2 = redlength;
2546 den = tstop + redlength;
2547 while ( tstop[len1-1] == 0 ) len1--;
2548 while ( den[len2-1] == 0 ) len2--;
2549 if ( GCDlen == 1 && GCDbuffer[0] == 1 ) {}
2551 GcdLong(BHEAD (UWORD *)tstop,len1,GCDbuffer,GCDlen,GCDbuffer2,&GCDlen2);
2552 ap = GCDbuffer; GCDbuffer = GCDbuffer2; GCDbuffer2 = ap;
2553 a = GCDlen; GCDlen = GCDlen2; GCDlen2 = a;
2555 if ( len2 == 1 && den[0] == 1 ) {}
2557 GcdLong(BHEAD LCMbuffer,LCMlen,(UWORD *)den,len2,LCMbuffer2,&LCMlen2);
2558 DivLong((UWORD *)den,len2,LCMbuffer2,LCMlen2,
2559 GCDbuffer2,&GCDlen2,(UWORD *)AT.WorkPointer,&a);
2560 MulLong(LCMbuffer,LCMlen,GCDbuffer2,GCDlen2,LCMbuffer2,&LCMlen2);
2561 ap = LCMbuffer; LCMbuffer = LCMbuffer2; LCMbuffer2 = ap;
2562 a = LCMlen; LCMlen = LCMlen2; LCMlen2 = a;
2566 if ( GCDlen != 1 || GCDbuffer[0] != 1 || LCMlen != 1 || LCMbuffer[0] != 1 ) {
2567 redlength = GCDlen;
if ( LCMlen > GCDlen ) redlength = LCMlen;
2568 for ( i = 0; i < GCDlen; i++ ) *tout++ = (WORD)(GCDbuffer[i]);
2569 for ( ; i < redlength; i++ ) *tout++ = 0;
2570 for ( i = 0; i < LCMlen; i++ ) *tout++ = (WORD)(LCMbuffer[i]);
2571 for ( ; i < redlength; i++ ) *tout++ = 0;
2572 *tout++ = (2*redlength+1)*sign;
2573 for ( i = 0; i < LCMlen; i++ ) *tt++ = (WORD)(LCMbuffer[i]);
2574 for ( ; i < redlength; i++ ) *tt++ = 0;
2575 for ( i = 0; i < GCDlen; i++ ) *tt++ = (WORD)(GCDbuffer[i]);
2576 for ( ; i < redlength; i++ ) *tt++ = 0;
2577 *tt++ = (2*redlength+1)*sign;
2581 *tout++ = 1; *tout++ = 1; *tout++ = 3*sign;
2582 *tt++ = 1; *tt++ = 1; *tt++ = 3*sign;
2583 if ( sign != 1 ) action++;
2585 NumberFree(LCMbuffer2,
"MakeInteger");
2586 NumberFree(LCMbuffer ,
"MakeInteger");
2587 NumberFree(GCDbuffer2,
"MakeInteger");
2588 NumberFree(GCDbuffer ,
"MakeInteger");
2595 *term = tout - term; *tout = 0;
2596 *tterm = tt - tterm; *tt = 0;
2597 AT.WorkPointer = tt;
2598 inp = MultiplyWithTerm(BHEAD in,tterm,2);
2599 AT.WorkPointer = tterm;
2600 t = inp;
while ( *t ) t += *t;
2601 j = (t-inp); t = inp;
2602 if ( j*
sizeof(WORD) > (size_t)(AM.MaxTer) )
goto OverWork;
2603 in = tout = TermMalloc(
"TakeSymbolContent");
2604 NCOPY(tout,t,j); *tout = 0;
2605 if ( AN.tryterm > 0 ) { TermFree(inp,
"MultiplyWithTerm"); AN.tryterm = 0; }
2606 else { M_free(inp,
"MultiplyWithTerm"); }
2609 t = in;
while ( *t ) t += *t;
2611 if ( j*
sizeof(WORD) > (size_t)(AM.MaxTer) )
goto OverWork;
2612 in = tout = TermMalloc(
"TakeSymbolContent");
2613 NCOPY(tout,t,j); *tout = 0;
2614 term[0] = 4; term[1] = 1; term[2] = 1; term[3] = 3; term[4] = 0;
2620 AT.WorkPointer = oldworkpointer;
2624 MLOCK(ErrorMessageLock);
2625 MesPrint(
"Term too complex. Maybe increasing MaxTermSize can help");
2626 MUNLOCK(ErrorMessageLock);
2628 MLOCK(ErrorMessageLock);
2629 MesCall(
"TakeSymbolContent");
2630 MUNLOCK(ErrorMessageLock);
2644void GCDclean(PHEAD WORD *num, WORD *den)
2646 WORD *out1 = TermMalloc(
"GCDclean");
2647 WORD *out2 = TermMalloc(
"GCDclean");
2648 WORD *t1, *t2, *r1, *r2, *t1stop, *t2stop, csize1, csize2, csize3, pow, sign;
2650 t1stop = num+*num; sign = ( t1stop[-1] < 0 ) ? -1 : 1;
2651 csize1 = ABS(t1stop[-1]); t1stop -= csize1;
2652 t2stop = den+*den;
if ( t2stop[-1] < 0 ) sign = -sign;
2653 csize2 = ABS(t2stop[-1]); t2stop -= csize2;
2654 t1 = num+1; t2 = den+1;
2655 r1 = out1+3; r2 = out2+3;
2656 if ( t1 == t1stop ) {
2657 if ( t2 < t2stop ) {
2658 for ( i = 2; i < t2[1]; i += 2 ) {
2659 if ( t2[i+1] < 0 ) { *r1++ = t2[i]; *r1++ = -t2[i+1]; }
2660 else { *r2++ = t2[i]; *r2++ = t2[i+1]; }
2664 else if ( t2 == t2stop ) {
2665 for ( i = 2; i < t1[1]; i += 2 ) {
2666 if ( t1[i+1] < 0 ) { *r2++ = t1[i]; *r2++ = -t1[i+1]; }
2667 else { *r1++ = t1[i]; *r1++ = t1[i+1]; }
2672 while ( t1 < t1stop && t2 < t2stop ) {
2674 if ( t1[1] > 0 ) { *r1++ = *t1; *r1++ = t1[1]; t1 += 2; }
2675 else if ( t1[1] < 0 ) { *r2++ = *t1; *r2++ = -t1[1]; t1 += 2; }
2677 else if ( *t1 > *t2 ) {
2678 if ( t2[1] > 0 ) { *r2++ = *t2; *r2++ = t2[1]; t2 += 2; }
2679 else if ( t2[1] < 0 ) { *r1++ = *t2; *r1++ = -t2[1]; t2 += 2; }
2683 if ( pow > 0 ) { *r1++ = *t1; *r1++ = pow; }
2684 else if ( pow < 0 ) { *r2++ = *t1; *r2++ = -pow; }
2688 while ( t1 < t1stop ) {
2689 if ( t1[1] < 0 ) { *r2++ = *t1; *r2++ = -t1[1]; }
2690 else { *r1++ = *t1; *r1++ = t1[1]; }
2693 while ( t2 < t2stop ) {
2694 if ( t2[1] < 0 ) { *r1++ = *t2; *r1++ = -t2[1]; }
2695 else { *r2++ = *t2; *r2++ = t2[1]; }
2699 if ( r1 > out1+3 ) { out1[1] = SYMBOL; out1[2] = r1 - out1 - 1; }
2701 if ( r2 > out2+3 ) { out2[1] = SYMBOL; out2[2] = r2 - out2 - 1; }
2706 csize1 = REDLENG(csize1);
2707 csize2 = REDLENG(csize2);
2708 if ( DivRat(BHEAD (UWORD *)t1stop,csize1,(UWORD *)t2stop,csize2,(UWORD *)r1,&csize3) ) {
2709 MLOCK(ErrorMessageLock);
2710 MesCall(
"GCDclean");
2711 MUNLOCK(ErrorMessageLock);
2714 UnPack((UWORD *)r1,csize3,&csize2,&csize1);
2715 t2 = r1+ABS(csize3);
2716 for ( i = 0; i < csize2; i++ ) r2[i] = t2[i];
2717 r2 += csize2; *r2++ = 1;
2718 for ( i = 1; i < csize2; i++ ) *r2++ = 0;
2719 csize2 = INCLENG(csize2); *r2++ = csize2; *out2 = r2-out2;
2720 r1 += ABS(csize1); *r1++ = 1;
2721 for ( i = 1; i < ABS(csize1); i++ ) *r1++ = 0;
2722 csize1 = INCLENG(csize1); *r1++ = csize1; *out1 = r1-out1;
2724 t1 = num; t2 = out1; i = *out1; NCOPY(t1,t2,i); *t1 = 0;
2725 if ( sign < 0 ) t1[-1] = -t1[-1];
2726 t1 = den; t2 = out2; i = *out2; NCOPY(t1,t2,i); *t1 = 0;
2728 TermFree(out2,
"GCDclean");
2729 TermFree(out1,
"GCDclean");
2741WORD *PolyDiv(PHEAD WORD *a,WORD *b,
char *text)
2747 return(poly_div(BHEAD a,b,1));
2786WORD divrem[4] = { DIVFUNCTION, REMFUNCTION, INVERSEFUNCTION, MULFUNCTION };
2788int DIVfunction(PHEAD WORD *term,WORD level,
int par)
2791 WORD *t, *tt, *r, *arg1 = 0, *arg2 = 0, *arg3 = 0, *termout;
2792 WORD *tstop, *tend, *r3, *rr, *rstop, tlength, rlength, newlength;
2793 WORD *proper1, *proper2, *proper3 = 0, numdol = -1;
2794 int numargs = 0, type1, type2, actionflag1, actionflag2;
2795 WORD startebuf = cbuf[AT.ebufnum].numrhs;
2796 int division = ( par <= 2 );
2797 if ( par < 0 || par > 3 ) {
2798 MLOCK(ErrorMessageLock);
2799 MesPrint(
"Internal error. Illegal parameter %d in DIVfunction.",par);
2800 MUNLOCK(ErrorMessageLock);
2806 tend = term + *term; tstop = tend - ABS(tend[-1]);
2808 while ( t < tstop ) {
2809 if ( *t != divrem[par] ) { t += t[1];
continue; }
2811 tt = t + t[1]; numargs = 0;
2813 if ( numargs == 0 ) { arg1 = r; }
2814 if ( numargs == 1 ) { arg2 = r; }
2815 if ( numargs == 2 && *r == -DOLLAREXPRESSION ) { numdol = r[1]; }
2819 if ( numargs == 2 )
break;
2820 if ( division && numargs == 3 )
break;
2824 MLOCK(ErrorMessageLock);
2825 MesPrint(
"Internal error. Indicated div_ or rem_ function not encountered.");
2826 MUNLOCK(ErrorMessageLock);
2832 if ( division && *arg1 == -SNUMBER && arg1[1] == 0 ) {
2833 if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2835 MLOCK(ErrorMessageLock);
2836 MesPrint(
"0/0 in either div_ or rem_ function.");
2837 MUNLOCK(ErrorMessageLock);
2840 if ( numdol >= 0 ) PutTermInDollar(0,numdol);
2843 if ( division && *arg2 == -SNUMBER && arg2[1] == 0 ) {
2845 MLOCK(ErrorMessageLock);
2846 MesPrint(
"Division by zero in either div_ or rem_ function.");
2847 MUNLOCK(ErrorMessageLock);
2851 if ( (*arg1 == -SNUMBER && arg1[1] == 0) ||
2852 (*arg2 == -SNUMBER && arg2[1] == 0) ) {
2856 if ( ( arg1 = ConvertArgument(BHEAD arg1, &type1) ) == 0 )
goto CalledFrom;
2857 if ( ( arg2 = ConvertArgument(BHEAD arg2, &type2) ) == 0 )
goto CalledFrom;
2858 if ( division && *arg1 == 0 ) {
2860 M_free(arg2,
"DIVfunction");
2861 M_free(arg1,
"DIVfunction");
2864 M_free(arg2,
"DIVfunction");
2865 M_free(arg1,
"DIVfunction");
2866 if ( numdol >= 0 ) PutTermInDollar(0,numdol);
2869 if ( division && *arg2 == 0 ) {
2870 M_free(arg2,
"DIVfunction");
2871 M_free(arg1,
"DIVfunction");
2874 if ( !division && (*arg1 == 0 || *arg2 == 0) ) {
2875 M_free(arg2,
"DIVfunction");
2876 M_free(arg1,
"DIVfunction");
2879 if ( ( proper1 = PutExtraSymbols(BHEAD arg1,startebuf,&actionflag1) ) == 0 )
goto CalledFrom;
2880 if ( ( proper2 = PutExtraSymbols(BHEAD arg2,startebuf,&actionflag2) ) == 0 )
goto CalledFrom;
2889 M_free(arg2,
"DIVfunction");
2898 M_free(arg1,
"DIVfunction");
2899 if ( par == 0 ) proper3 = poly_div(BHEAD proper1, proper2,0);
2900 else if ( par == 1 ) proper3 = poly_rem(BHEAD proper1, proper2,0);
2901 else if ( par == 2 ) proper3 = poly_inverse(BHEAD proper1, proper2);
2902 else if ( par == 3 ) proper3 = poly_mul(BHEAD proper1, proper2);
2903 if ( proper3 == 0 )
goto CalledFrom;
2904 if ( actionflag1 || actionflag2 ) {
2905 if ( ( arg3 = TakeExtraSymbols(BHEAD proper3,startebuf) ) == 0 )
goto CalledFrom;
2906 M_free(proper3,
"DIVfunction");
2911 M_free(proper2,
"DIVfunction");
2912 M_free(proper1,
"DIVfunction");
2913 cbuf[AT.ebufnum].numrhs = startebuf;
2915 termout = AT.WorkPointer;
2917 tlength = REDLENG(tlength);
2920 tt = term + 1; rr = termout + 1;
2921 while ( tt < t ) *rr++ = *tt++;
2924 rstop = r3 - ABS(r3[-1]);
2925 while ( r < rstop ) *rr++ = *r++;
2927 while ( tt < tstop ) *rr++ = *tt++;
2929 rlength = REDLENG(rlength);
2930 if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)rstop,rlength,
2931 (UWORD *)rr,&newlength) < 0 )
goto CalledFrom;
2932 rlength = INCLENG(newlength);
2935 *termout = rr - termout;
2936 AT.WorkPointer = rr;
2937 if (
Generator(BHEAD termout,level) )
goto CalledFrom;
2939 AT.WorkPointer = termout;
2941 M_free(arg3,
"DIVfunction");
2944 MLOCK(ErrorMessageLock);
2945 MesCall(
"DIVfunction");
2946 MUNLOCK(ErrorMessageLock);
2957WORD *MULfunc(PHEAD WORD *p1, WORD *p2)
2959 WORD *prod,size1,size2,size3,*t,*tfill,*ps1,*ps2,sign1,sign2, error, *p3;
2960 UWORD *num1, *num2, *num3;
2962 WORD oldsorttype = AR.SortType;
2963 COMPARE oldcompareroutine = AR.CompareRoutine;
2964 AR.SortType = SORTHIGHFIRST;
2966 num3 = NumberMalloc(
"MULfunc");
2967 prod = TermMalloc(
"MULfunc");
2970 ps1 = p1+*p1; num1 = (UWORD *)(ps1 - ABS(ps1[-1])); size1 = ps1[-1];
2971 if ( size1 < 0 ) { sign1 = -1; size1 = -size1; }
2973 size1 = (size1-1)/2;
2976 ps2 = p3+*p3; num2 = (UWORD *)(ps2 - ABS(ps2[-1])); size2 = ps2[-1];
2977 if ( size2 < 0 ) { sign2 = -1; size2 = -size2; }
2979 size2 = (size2-1)/2;
2980 if ( MulLong(num1,size1,num2,size2,num3,&size3) ) {
2983 MLOCK(ErrorMessageLock);
2984 MesPrint(
" Error %d",error);
2986 MUNLOCK(ErrorMessageLock);
2990 t = p1+1;
while ( t < (WORD *)num1 ) *tfill++ = *t++;
2991 t = p3+1;
while ( t < (WORD *)num2 ) *tfill++ = *t++;
2993 for ( i = 0; i < size3; i++ ) *tfill++ = *t++;
2995 for ( i = 1; i < size3; i++ ) *tfill++ = 0;
2996 *tfill++ = (2*size3+1)*sign1*sign2;
2997 prod[0] = tfill - prod;
2999 if (
StoreTerm(BHEAD prod) ) { error = 3;
goto CalledFrom; }
3004 NumberFree(num3,
"MULfunc");
3006 AR.CompareRoutine = oldcompareroutine;
3007 AR.SortType = oldsorttype;
3018WORD *ConvertArgument(PHEAD WORD *arg,
int *type)
3020 WORD *output, *t, *r;
3023 output = (WORD *)Malloc1((*arg)*
sizeof(WORD),
"ConvertArgument");
3024 i = *arg - ARGHEAD; t = arg + ARGHEAD; r = output;
3029 if ( *arg == -EXPRESSION ) {
3031 return(CreateExpression(BHEAD arg[1]));
3033 if ( *arg == -DOLLAREXPRESSION ) {
3036 d = DolToTerms(BHEAD arg[1]);
3042 output = (WORD *)Malloc1((d->size+1)*
sizeof(WORD),
"Copy of dollar content");
3043 WCOPY(output,d->where,d->size+1);
3044 if ( d->factors ) { M_free(d->factors,
"Dollar factors"); d->factors = 0; }
3045 M_free(d,
"Copy of dollar variable");
3053 output = (WORD *)Malloc1(size*
sizeof(WORD),
"ConvertArgument");
3056 output[0] = 8; output[1] = SYMBOL; output[2] = 4; output[3] = arg[1];
3057 output[4] = 1; output[5] = 1; output[6] = 1; output[7] = 3;
3063 output[0] = 7; output[1] = INDEX; output[2] = 3; output[3] = arg[1];
3064 output[4] = 1; output[5] = 1;
3065 if ( *arg == -MINVECTOR ) output[6] = -3;
3072 output[1] = -arg[1]; output[2] = 1; output[3] = -3;
3075 output[1] = arg[1]; output[2] = 1; output[3] = 3;
3080 output[0] = FUNHEAD+4;
3082 output[2] = FUNHEAD;
3083 for ( i = 3; i <= FUNHEAD; i++ ) output[i] = 0;
3084 output[FUNHEAD+1] = 1;
3085 output[FUNHEAD+2] = 1;
3086 output[FUNHEAD+3] = 3;
3087 output[FUNHEAD+4] = 0;
3105char *TheErrorMessage[] = {
3106 "PolyRatFun not of a type that FORM will expand: incorrect variable inside."
3107 ,
"Division by zero in PolyRatFun encountered in ExpandRat."
3108 ,
"Irregular code in PolyRatFun encountered in ExpandRat."
3109 ,
"Called from ExpandRat."
3110 ,
"WorkSpace overflow. Change parameter WorkSpace in setup file?"
3113int ExpandRat(PHEAD WORD *fun)
3115 WORD *r, *rr, *rrr, *tt, *tnext, *arg1, *arg2, *rmin = 0, *rmininv;
3116 WORD *rcoef, rsize, rcopy, *ow = AT.WorkPointer;
3117 WORD *numerator, *denominator, *rnext;
3118 WORD *thecopy, *rc, ncoef, newcoef, *m, *mm, nco, *outarg = 0;
3119 UWORD co[2], co1[2], co2[2];
3120 WORD OldPolyFunPow = AR.PolyFunPow;
3121 int i, j, minpow = 0, eppow, first, error = 0, ipoly;
3122 if ( fun[1] == FUNHEAD ) {
return(0); }
3123 tnext = fun + fun[1];
3124 if ( fun[1] == fun[FUNHEAD]+FUNHEAD ) {
3125 if ( fun[2] == 0 ) {
goto done; }
3130 if ( outarg == 0 ) outarg = TermMalloc(
"ExpandRat")+ARGHEAD;
3133 r = fun+FUNHEAD+ARGHEAD;
3134 if ( AR.PolyFunExp == 2 ) {
3135 WORD minpow2 = MAXPOWER, *rrm;
3137 while ( rrm < tnext ) {
3139 if ( minpow2 > 0 ) minpow2 = 0;
3141 else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3142 if ( minpow2 > 0 ) minpow2 = 0;
3145 if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3146 if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3149 MesPrint(
"Illegal term in expanded polyratfun.");
3155 AR.PolyFunPow += minpow2;
3157 while ( r < tnext ) {
3159 i = *r; rrr = outarg; NCOPY(rrr,r,i);
3160 Normalize(BHEAD outarg);
3161 if ( *outarg > 0 )
StoreTerm(BHEAD outarg);
3163 r = fun+FUNHEAD+ARGHEAD;
3167 fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0;
3172 if ( ToFast(rr,rr) ) {
3173 NEXTARG(rr); fun[1] = rr - fun;
3176 while ( *r ) r += *r;
3177 *rr = r-rr; rr[1] = CLEANFLAG;
3190 arg1 = tt; NEXTARG(tt);
3192 arg2 = tt; NEXTARG(tt);
3193 if ( tt != tnext ) { arg1 = arg2 = 0; }
3197 if ( *arg1 < 0 ) { fun[2] = CLEANFLAG;
goto done; }
3198 if ( fun[2] == CLEANFLAG )
goto done;
3204 if ( outarg == 0 ) outarg = TermMalloc(
"ExpandRat")+ARGHEAD;
3213 if ( *arg2 == -SYMBOL && arg2[1] == AR.PolyFunVar ) {
3214 rr = r = fun+FUNHEAD+ARGHEAD;
3216 if ( *arg1 == -SYMBOL ) {
3217 if ( arg1[1] == AR.PolyFunVar ) {
3218 *r++ = 4; *r++ = 1; *r++ = 1; *r++ = 3; *r = 0;
3221 *r++ = 10; *r++ = SYMBOL; *r++ = 6;
3222 *r++ = arg1[1]; *r++ = 1;
3223 *r++ = AR.PolyFunVar; *r++ = -1;
3224 *r++ = 1; *r++ = 1; *r++ = 3; *r = 0;
3225 Normalize(BHEAD rr);
3228 else if ( *arg1 == -SNUMBER ) {
3230 if ( nco == 0 ) { *r++ = 0; }
3232 *r++ = 8; *r++ = SYMBOL; *r++ = 4;
3233 *r++ = AR.PolyFunVar; *r++ = -1;
3234 *r++ = ABS(nco); *r++ = 1;
3235 if ( nco < 0 ) *r++ = -3;
3240 else { error = 2;
goto onerror; }
3245 if ( AR.PolyFunExp == 2 ) {
3246 WORD minpow2 = MAXPOWER, *rrm;
3248 while ( rrm < arg2 ) {
3250 if ( minpow2 > 0 ) minpow2 = 0;
3252 else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3253 if ( minpow2 > 0 ) minpow2 = 0;
3256 if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3257 if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3260 MesPrint(
"Illegal term in expanded polyratfun.");
3266 AR.PolyFunPow += minpow2-1;
3268 while ( m < arg2 ) {
3270 rrr = r++; mm = m + *m;
3271 *r++ = SYMBOL; *r++ = 4; *r++ = AR.PolyFunVar; *r++ = -1;
3272 m++;
while ( m < mm ) *r++ = *m++;
3274 Normalize(BHEAD rrr);
3278 r = rr;
while ( *r ) r += *r;
3281 fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = CLEANFLAG;
3288 if ( ToFast(rr,rr) ) {
3292 else { fun[1] = r - fun; }
3297 else if ( *arg2 == -SNUMBER ) {
3299 if ( arg2[1] == 0 ) { error = 1;
goto onerror; }
3300 if ( *arg1 == -SNUMBER ) {
3301 if ( arg1[1] == 0 ) { *r++ = 0; }
3303 co1[0] = ABS(arg1[1]); co1[1] = 1;
3304 co2[0] = 1; co2[1] = ABS(arg2[1]);
3305 MulRat(BHEAD co1,1,co2,1,co,&nco);
3306 *r++ = 4; *r++ = (WORD)(co[0]); *r++ = (WORD)(co[1]);
3307 if ( ( arg1[1] < 0 && arg2[1] > 0 ) ||
3308 ( arg1[1] > 0 && arg2[1] < 0 ) ) *r++ = -3;
3313 else if ( *arg1 == -SYMBOL ) {
3314 *r++ = 8; *r++ = SYMBOL; *r++ = 4;
3315 *r++ = arg1[1]; *r++ = 1;
3316 *r++ = 1; *r++ = ABS(arg2[1]);
3317 if ( arg2[1] < 0 ) *r++ = -3;
3321 else if ( *arg1 < 0 ) { error = 2;
goto onerror; }
3325 if ( AR.PolyFunExp == 2 ) {
3326 WORD minpow2 = MAXPOWER, *rrm;
3328 while ( rrm < arg2 ) {
3330 if ( minpow2 > 0 ) minpow2 = 0;
3332 else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3333 if ( minpow2 > 0 ) minpow2 = 0;
3336 if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3337 if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3340 MesPrint(
"Illegal term in expanded polyratfun.");
3346 AR.PolyFunPow += minpow2;
3348 while ( m < arg2 ) {
3350 rrr = r++; mm = m + *m;
3351 *r++ = DENOMINATOR; *r++ = FUNHEAD + 2; *r++ = DIRTYFLAG;
3353 *r++ = arg2[0]; *r++ = arg2[1];
3354 m++;
while ( m < mm ) *r++ = *m++;
3356 if ( r < AT.WorkTop && r >= AT.WorkSpace )
3358 Normalize(BHEAD rrr);
3359 if ( ABS(rrr[*rrr-1]) == *rrr-1 ) {
3360 if ( AR.PolyFunPow >= 0 ) {
3364 else if ( rrr[1] == SYMBOL && rrr[2] == 4 &&
3365 rrr[3] == AR.PolyFunVar && rrr[4] <= AR.PolyFunPow ) {
3371 r = rr;
while ( *r ) r += *r;
3373 r = fun + FUNHEAD + ARGHEAD;
3376 *rr = r - rr; rr[1] = CLEANFLAG;
3377 if ( ToFast(rr,rr) ) {
3381 else { fun[1] = r - fun; }
3385 else { error = 0;
goto onerror; }
3390 while ( r < tnext ) {
3391 rr = r + *r; rr -= ABS(rr[-1]);
3393 if ( first ) { minpow = 0; first = 0; rmin = r; }
3394 else if ( minpow > 0 ) { minpow = 0; rmin = r; }
3396 else if ( r[1] != SYMBOL || r[2] != 4 || r[3] != AR.PolyFunVar
3397 || r[4] > MAXPOWER ) { error = 0;
goto onerror; }
3398 else if ( first ) { minpow = r[4]; first = 0; rmin = r; }
3399 else if ( r[4] < minpow ) { minpow = r[4]; rmin = r; }
3417 AT.WorkPointer += AM.MaxTer/
sizeof(WORD);
3418 if ( AT.WorkPointer + (AM.MaxTer/
sizeof(WORD)) >= AT.WorkTop ) {
3419 error = 4;
goto onerror;
3421 rmininv = r = AT.WorkPointer;
3422 rr = rmin; i = *rmin; NCOPY(r,rr,i)
3423 if ( minpow != 0 ) { rmininv[4] = -rmininv[4]; }
3426 rsize = (rsize-1)/2; rr = rcoef + rsize;
3427 for ( i = 0; i < rsize; i++ ) {
3428 rcopy = rcoef[i]; rcoef[i] = rr[i]; rr[i] = rcopy;
3432 ToGeneral(arg1,r,0);
3433 arg1 = r; r += *r; *r++ = 0; rcopy = 0;
3438 rcopy = *r; *r++ = 0;
3443 AT.LeaveNegative = 1;
3444 numerator = MultiplyWithTerm(BHEAD arg1+ARGHEAD,rmininv,0);
3445 AT.LeaveNegative = 0;
3447 r = numerator;
while ( *r ) r += *r;
3448 AT.WorkPointer = r+1;
3449 rcopy = arg2[*arg2]; arg2[*arg2] = 0;
3450 denominator = MultiplyWithTerm(BHEAD arg2+ARGHEAD,rmininv,1);
3451 arg2[*arg2] = rcopy;
3452 r = denominator;
while ( *r ) r += *r;
3453 AT.WorkPointer = r+1;
3460 rr = r + *r; rr -= ABS(rr[-1]);
3462 if ( first ) { minpow = 0; first = 0; }
3463 else if ( minpow > 0 ) { minpow = 0; }
3465 else if ( r[1] != SYMBOL ) { error = 0;
goto onerror; }
3467 for ( i = 3; i < r[2]; i += 2 ) {
3468 if ( r[i] == AR.PolyFunVar ) {
3469 if ( first ) { minpow = r[i+1]; first = 0; }
3470 else if ( r[i+1] < minpow ) minpow = r[i+1];
3475 if ( first ) { minpow = 0; first = 0; }
3476 else if ( minpow > 0 ) minpow = 0;
3486 if ( AR.PolyFunExp == 3 ) {
3487 ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow-minpow,AR.PolyFunVar);
3490 ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow,AR.PolyFunVar);
3502 rrr = rnext - ABS(rnext[-1]);
3506 j = rr[1] - 2; rr += 2;
3508 if ( *rr == AR.PolyFunVar ) { eppow = rr[1];
break; }
3515 for ( i = 0; i <= AR.PolyFunPow-eppow+minpow; i++ ) {
3516 if ( AT.pWorkSpace[ipoly+i] == 0 )
continue;
3521 rr = thecopy = AT.WorkPointer;
3522 while ( rc < rrr ) *rr++ = *rc++;
3524 *rr++ = SYMBOL; *rr++ = 4; *rr++ = AR.PolyFunVar; *rr++ = i;
3526 ncoef = REDLENG(rnext[-1]);
3527 MulRat(BHEAD (UWORD *)rrr,ncoef,
3528 (UWORD *)(AT.pWorkSpace[ipoly+i])+1,AT.pWorkSpace[ipoly+i][0]
3529 ,(UWORD *)rr,&newcoef);
3530 ncoef = ABS(newcoef); rr += 2*ncoef;
3531 newcoef = INCLENG(newcoef);
3533 *thecopy = rr - thecopy;
3534 AT.WorkPointer = rr;
3535 Normalize(BHEAD thecopy);
3536 if ( *thecopy > 0 )
StoreTerm(BHEAD thecopy);
3537 AT.WorkPointer = thecopy;
3544 rr = fun + FUNHEAD; r = rr + ARGHEAD;
3547 fun[1] = FUNHEAD+2; fun[2] = CLEANFLAG;
3548 fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0;
3551 while ( *r ) r += *r;
3552 rr[0] = r-rr; rr[1] = CLEANFLAG;
3553 if ( ToFast(rr,rr) ) { NEXTARG(rr); fun[1] = rr-fun; }
3554 else { fun[1] = r-fun; }
3559 if ( outarg ) TermFree(outarg-ARGHEAD,
"ExpandRat");
3560 AR.PolyFunPow = OldPolyFunPow;
3561 AT.WorkPointer = ow;
3562 AN.PolyNormFlag = 1;
3565 if ( outarg ) TermFree(outarg-ARGHEAD,
"ExpandRat");
3566 AR.PolyFunPow = OldPolyFunPow;
3567 AT.WorkPointer = ow;
3568 MLOCK(ErrorMessageLock);
3569 MesPrint(TheErrorMessage[error]);
3570 MUNLOCK(ErrorMessageLock);
3592int InvPoly(PHEAD WORD *inpoly, WORD maxpow, WORD sym)
3594 int needed, inpointers, outpointers, maxinput = 0, i, j;
3595 WORD *t, *tt, *ttt, *w, *c, *cc, *ccc, lenc, lenc1, lenc2, rc, *c1, *c2;
3599 needed = (maxpow+1)*2;
3600 WantAddPointers(needed);
3601 inpointers = AT.pWorkPointer;
3602 outpointers = AT.pWorkPointer+maxpow+1;
3603 for ( i = 0; i < needed; i++ ) AT.pWorkSpace[inpointers+i] = 0;
3613 if ( t[1] != 1 || t[2] != 1 || t[3] != 3 )
goto onerror;
3614 AT.pWorkSpace[inpointers] = 0;
3616 else if ( t[1] != SYMBOL || t[2] != 4 || t[3] != sym || t[4] < 0 )
goto onerror;
3617 else if ( t[4] > maxpow ) {}
3619 if ( t[4] > maxinput ) maxinput = t[4];
3620 AT.pWorkSpace[inpointers+t[4]] = w;
3621 tt = t + *t; rc = -*--tt;
3622 rc = REDLENG(rc); *w++ = rc;
3624 while ( ttt < tt ) *w++ = *ttt++;
3632 AT.pWorkSpace[outpointers] = w;
3633 *w++ = 1; *w++ = 1; *w++ = 1;
3634 c = TermMalloc(
"InvPoly");
3635 c1 = TermMalloc(
"InvPoly");
3636 c2 = TermMalloc(
"InvPoly");
3637 for ( j = 1; j <= maxpow; j++ ) {
3641 if ( ( cc = AT.pWorkSpace[inpointers+j] ) != 0 ) {
3643 i = 2*ABS(lenc); ccc = c;
3647 for ( i = MiN(j-1,maxinput); i > 0; i-- ) {
3651 if ( AT.pWorkSpace[inpointers+i] == 0
3652 || AT.pWorkSpace[outpointers+j-i] == 0 ) {
3655 if ( MulRat(BHEAD (UWORD *)(AT.pWorkSpace[inpointers+i]+1),AT.pWorkSpace[inpointers+i][0],
3656 (UWORD *)(AT.pWorkSpace[outpointers+j-i]+1),AT.pWorkSpace[outpointers+j-i][0],
3657 (UWORD *)c1,&lenc1) )
goto calcerror;
3659 cc = c; c = c1; c1 = cc;
3663 if ( AddRat(BHEAD (UWORD *)c,lenc,(UWORD *)c1,lenc1,(UWORD *)c2,&lenc2) )
3665 cc = c; c = c2; c2 = cc;
3673 if ( lenc == 0 ) AT.pWorkSpace[outpointers+j] = 0;
3675 AT.pWorkSpace[outpointers+j] = w;
3677 i = 2*ABS(lenc); ccc = c;
3682 TermFree(c2,
"InvPoly");
3683 TermFree(c1,
"InvPoly");
3684 TermFree(c ,
"InvPoly");
3686 return(outpointers);
3688 MLOCK(ErrorMessageLock);
3689 MesPrint(
"Incorrect symbol field in InvPoly.");
3690 MUNLOCK(ErrorMessageLock);
3694 MLOCK(ErrorMessageLock);
3695 MesPrint(
"Called from InvPoly.");
3696 MUNLOCK(ErrorMessageLock);
WORD CompareSymbols(WORD *, WORD *, WORD)
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
LONG EndSort(PHEAD WORD *, int)
WORD InsertTerm(PHEAD WORD *, WORD, WORD, WORD *, WORD *, WORD)
WORD Generator(PHEAD WORD *, WORD)
WORD StoreTerm(PHEAD WORD *)
int GetModInverses(WORD, WORD, WORD *, WORD *)
int SymbolNormalize(WORD *)
WORD * poly_gcd(PHEAD WORD *, WORD *, WORD)
WORD * TakeContent(PHEAD WORD *in, WORD *term)
WORD * TakeSymbolContent(PHEAD WORD *in, WORD *term)