62 WORD 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];
93 We 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 );
170 WORD 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);
320 WORD 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;
395 WORD 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);
461 WORD 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");
520 WORD 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);
609 int 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, *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]; }
680 else if ( *tf != firstshort ) {
681 if ( *tf != -INDEX && *tf != -VECTOR && *t != -MINVECTOR ) {
684 if ( firstshort != -INDEX && firstshort != -VECTOR && firstshort != -MINVECTOR ) {
687 if ( tf[1] != firstvalue ) {
690 if ( *t == -MINVECTOR ) { firstshort = -VECTOR; }
691 if ( firstshort == -MINVECTOR ) { firstshort = -VECTOR; }
693 else if ( *tf > -FUNCTION && *tf != -SNUMBER && tf[1] != firstvalue ) {
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;
710 if ( firstvalue == 1 ) { gcdisone = 1;
break; }
714 termout = AT.WorkPointer;
715 AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
716 if ( AT.WorkPointer > AT.WorkTop ) {
717 MLOCK(ErrorMessageLock);
719 MUNLOCK(ErrorMessageLock);
728 i = t - term; tin = term; tout = termout;
730 if ( gcdisone || ( firstshort == -SNUMBER && firstvalue == 1 ) ) {
733 tin += t[1]; tstop = term + *term;
734 while ( tin < tstop ) *tout++ = *tin++;
735 *termout = tout - termout;
736 if ( sign < 0 ) tout[-1] = -tout[-1];
737 AT.WorkPointer = tout;
738 if (
Generator(BHEAD termout,level) < 0 )
goto CalledFrom;
739 AT.WorkPointer = termout;
740 AT.pWorkPointer = args;
747 if ( numargs == 0 ) {
750 if ( firstshort == 0 )
goto gcdone;
751 if ( firstshort == -SNUMBER ) {
752 *tout++ = SNUMBER; *tout++ = 4; *tout++ = firstvalue; *tout++ = 1;
755 else if ( firstshort == -SYMBOL ) {
756 *tout++ = SYMBOL; *tout++ = 4; *tout++ = firstvalue; *tout++ = 1;
759 else if ( firstshort == -VECTOR || firstshort == -INDEX ) {
760 *tout++ = INDEX; *tout++ = 3; *tout++ = firstvalue;
goto gcdone;
762 else if ( firstshort == -MINVECTOR ) {
764 *tout++ = INDEX; *tout++ = 3; *tout++ = firstvalue;
goto gcdone;
766 else if ( firstshort <= -FUNCTION ) {
767 *tout++ = firstvalue; *tout++ = FUNHEAD; FILLFUN(tout);
771 MLOCK(ErrorMessageLock);
772 MesPrint(
"Internal error. Illegal short argument in GCDfunction.");
773 MUNLOCK(ErrorMessageLock);
785 switch ( firstshort ) {
787 sh[0] = 4; sh[1] = ABS(firstvalue); sh[2] = 1;
788 if ( firstvalue < 0 ) sh[3] = -3;
795 sh[0] = 8; sh[1] = INDEX; sh[2] = 3; sh[3] = firstvalue;
796 sh[4] = 1; sh[5] = 1;
797 if ( firstshort == -MINVECTOR ) sh[6] = -3;
802 sh[0] = 8; sh[1] = SYMBOL; sh[2] = 4; sh[3] = firstvalue; sh[4] = 1;
803 sh[5] = 1; sh[6] = 1; sh[7] = 3; sh[8] = 0;
806 sh[0] = FUNHEAD+4; sh[1] = firstshort; sh[2] = FUNHEAD;
807 for ( i = 2; i < FUNHEAD; i++ ) sh[i+1] = 0;
808 sh[FUNHEAD+1] = 1; sh[FUNHEAD+2] = 1; sh[FUNHEAD+3] = 3; sh[FUNHEAD+4] = 0;
819 for ( i = 1; i < numargs; i++ ) {
820 for ( ii = i; ii > 0; ii-- ) {
821 arg1 = AT.pWorkSpace[args+ii];
822 arg2 = AT.pWorkSpace[args+ii-1];
825 if ( *arg1 == -EXPRESSION )
break;
826 if ( *arg2 == -DOLLAREXPRESSION )
break;
827 AT.pWorkSpace[args+ii] = arg2;
828 AT.pWorkSpace[args+ii-1] = arg1;
832 else if ( *arg2 < 0 ) {
833 AT.pWorkSpace[args+ii] = arg2;
834 AT.pWorkSpace[args+ii-1] = arg1;
837 if ( *arg1 > *arg2 ) {
838 AT.pWorkSpace[args+ii] = arg2;
839 AT.pWorkSpace[args+ii-1] = arg1;
852 for ( i = istart; i < numargs; i++ ) {
853 arg1 = AT.pWorkSpace[args+i];
855 oldval1 = arg1[*arg1]; arg1[*arg1] = 0;
858 GCDterms(BHEAD mh,m,mh); m += *m;
859 if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
860 gcdisone = 1; sign = 1; arg1[*arg1] = oldval1;
goto gcdone;
863 arg1[*arg1] = oldval1;
865 else if ( *arg1 == -DOLLAREXPRESSION ) {
866 if ( ( d = DolToTerms(BHEAD arg1[1]) ) != 0 ) {
869 GCDterms(BHEAD mh,m,mh); m += *m;
870 if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
871 gcdisone = 1; sign = 1;
872 if ( d->factors ) M_free(d->factors,
"Dollar factors");
873 M_free(d,
"Copy of dollar variable");
goto gcdone;
876 if ( d->factors ) M_free(d->factors,
"Dollar factors");
877 M_free(d,
"Copy of dollar variable");
881 mm = CreateExpression(BHEAD arg1[1]);
884 GCDterms(BHEAD mh,m,mh); m += *m;
885 if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
886 gcdisone = 1; sign = 1; M_free(mm,
"CreateExpression");
goto gcdone;
889 M_free(mm,
"CreateExpression");
894 firstshort = -SNUMBER; firstvalue = mh[1] * (mh[3]/3);
896 else if ( mh[1] == SYMBOL ) {
897 firstshort = -SYMBOL; firstvalue = mh[3];
899 else if ( mh[1] == INDEX ) {
900 firstshort = -INDEX; firstvalue = mh[3];
901 if ( mh[6] == -3 ) firstshort = -MINVECTOR;
903 else if ( mh[1] >= FUNCTION ) {
904 firstshort = -mh[1]; firstvalue = mh[1];
924 for ( i = 0; i < numargs; i++ ) {
925 arg1 = AT.pWorkSpace[args+i];
926 if ( *arg1 > 0 && arg1[ARGHEAD]+ARGHEAD == *arg1 ) {
931 arg2 = AT.pWorkSpace[args];
932 AT.pWorkSpace[args] = arg1;
933 AT.pWorkSpace[args+1] = arg2;
935 m = mh = AT.WorkPointer;
936 mm = arg1+ARGHEAD; i = *mm;
956 for ( i = 0; i < numargs; i++ ) {
957 arg1 = AT.pWorkSpace[args+i];
959 m = (WORD *)Malloc1(*arg1*
sizeof(WORD),
"argbuffer type 0");
968 else if ( *arg1 == -DOLLAREXPRESSION ) {
969 d = DolToTerms(BHEAD arg1[1]);
970 abuf[i].buffer = d->where;
973 m = abuf[i].buffer;
while ( *m ) m+= *m;
974 abuf[i].size = m-abuf[i].buffer;
976 else if ( *arg1 == -EXPRESSION ) {
977 abuf[i].buffer = CreateExpression(BHEAD arg1[1]);
979 m = abuf[i].buffer;
while ( *m ) m+= *m;
980 abuf[i].size = m-abuf[i].buffer;
983 MLOCK(ErrorMessageLock);
984 MesPrint(
"What argument is this?");
985 MUNLOCK(ErrorMessageLock);
989 for ( i = 0; i < numargs; i++ ) {
990 arg1 = abuf[i].buffer;
991 if ( arg1[*arg1] == 0 ) {
995 ab = abuf[i]; abuf[i] = abuf[0]; abuf[0] = ab;
997 for ( j = 1; j < numargs; j++ ) {
1000 GCDterms(BHEAD mh,m,mh); m += *m;
1001 if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
1002 gcdisone = 1; sign = 1;
break;
1007 mm = mh + *mh;
if ( mm[-1] < 0 ) { sign = -1; mm[-1] = -mm[-1]; }
1008 mstop = mm - mm[-1]; m = mh+1; mlength = mm[-1];
1009 while ( tin < t ) *tout++ = *tin++;
1010 while ( m < mstop ) *tout++ = *m++;
1012 while ( tin < tstop ) *tout++ = *tin++;
1013 tlength = REDLENG(tlength);
1014 mlength = REDLENG(mlength);
1015 if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)mstop,mlength,
1016 (UWORD *)tout,&newlength) < 0 )
goto CalledFrom;
1017 mlength = INCLENG(newlength);
1018 tout += ABS(mlength);
1019 tout[-1] = mlength*sign;
1020 *termout = tout - termout;
1021 AT.WorkPointer = tout;
1022 if (
Generator(BHEAD termout,level) < 0 )
goto CalledFrom;
1030 for ( i = 1; i < numargs; i++ ) {
1031 for ( ii = i; ii > 0; ii-- ) {
1032 if ( abuf[ii-1].size <= abuf[ii].size )
break;
1033 ab = abuf[ii-1]; abuf[ii-1] = abuf[ii]; abuf[ii] = ab;
1040 gcdout = abuf[0].buffer;
1041 for ( i = 1; i < numargs; i++ ) {
1042 g = GCDfunction3(BHEAD gcdout,abuf[i].buffer);
1043 if ( gcdout != abuf[0].buffer ) M_free(gcdout,
"gcdout");
1045 if ( gcdout[*gcdout] == 0 && gcdout[0] == 4 && gcdout[1] == 1
1046 && gcdout[2] == 1 && gcdout[3] == 3 )
break;
1050 tlength = REDLENG(tlength);
1052 tin = term; tout = termout;
while ( tin < t ) *tout++ = *tin++;
1054 mnext = mm + *mm; mlength = mnext[-1]; mstop = mnext - ABS(mlength);
1056 while ( mm < mstop ) *tout++ = *mm++;
1057 while ( tin < tstop ) *tout++ = *tin++;
1058 mlength = REDLENG(mlength);
1059 if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)mm,mlength,
1060 (UWORD *)tout,&newlength) < 0 )
goto CalledFrom;
1061 mlength = INCLENG(newlength);
1062 tout += ABS(mlength);
1064 *termout = tout - termout;
1065 AT.WorkPointer = tout;
1066 if (
Generator(BHEAD termout,level) < 0 )
goto CalledFrom;
1069 if ( action && ( gcdout != abuf[0].buffer ) ) M_free(gcdout,
"gcdout");
1076 for ( i = 0; i < numargs; i++ ) {
1077 if ( abuf[i].type == 0 ) { M_free(abuf[i].buffer,
"argbuffer type 0"); }
1078 else if ( abuf[i].type == 1 ) {
1080 if ( d->factors ) M_free(d->factors,
"Dollar factors");
1081 M_free(d,
"Copy of dollar variable");
1083 else if ( abuf[i].type == 2 ) { M_free(abuf[i].buffer,
"CreateExpression"); }
1085 M_free(abuf,
"argbuffer");
1090 AT.pWorkPointer = args;
1091 AT.WorkPointer = termout;
1095 MLOCK(ErrorMessageLock);
1096 MesCall(
"GCDfunction");
1097 MUNLOCK(ErrorMessageLock);
1122 WORD *GCDfunction3(
PHEAD WORD *in1, WORD *in2)
1125 WORD oldsorttype = AR.SortType;
1126 WORD *t, *tt, *gcdout, *term1, *term2, *confree1, *confree2, *gcdout1, *proper1, *proper2;
1127 int i, actionflag1, actionflag2;
1128 WORD startebuf = cbuf[AT.ebufnum].numrhs;
1129 if ( in2[*in2] == 0 ) { t = in1; in1 = in2; in2 = t; }
1130 if ( in1[*in1] == 0 ) {
1131 gcdout = (WORD *)Malloc1((*in1+1)*
sizeof(WORD),
"gcdout");
1132 i = *in1; t = gcdout; tt = in1; NCOPY(t,tt,i); *t = 0;
1135 GCDterms(BHEAD gcdout,t,gcdout);
1136 if ( gcdout[0] == 4 && gcdout[1] == 1
1137 && gcdout[2] == 1 && gcdout[3] == 3 )
break;
1146 AR.SortType = SORTHIGHFIRST;
1147 term1 = TermMalloc(
"GCDfunction3-a");
1148 term2 = TermMalloc(
"GCDfunction3-b");
1151 GCDterms(BHEAD term1,term2,term1);
1152 TermFree(term2,
"GCDfunction3-b");
1157 if ( ( proper1 = PutExtraSymbols(BHEAD confree1,startebuf,&actionflag1) ) == 0 )
goto CalledFrom;
1158 if ( confree1 != in1 ) M_free(confree1,
"TakeContent");
1159 if ( ( proper2 = PutExtraSymbols(BHEAD confree2,startebuf,&actionflag2) ) == 0 )
goto CalledFrom;
1160 if ( confree2 != in2 ) M_free(confree2,
"TakeContent");
1164 gcdout1 = poly_gcd(BHEAD proper1,proper2);
1165 M_free(proper1,
"PutExtraSymbols");
1166 M_free(proper2,
"PutExtraSymbols");
1168 AR.SortType = oldsorttype;
1169 if ( actionflag1 || actionflag2 ) {
1170 if ( ( gcdout = TakeExtraSymbols(BHEAD gcdout1,startebuf) ) == 0 )
goto CalledFrom;
1171 M_free(gcdout1,
"gcdout");
1177 cbuf[AT.ebufnum].numrhs = startebuf;
1181 if ( term1[0] != 4 || term1[3] != 3 || term1[1] != 1 || term1[2] != 1 ) {
1182 if ( ( gcdout1 = MultiplyWithTerm(BHEAD gcdout,term1) ) == 0 )
goto CalledFrom;
1183 M_free(gcdout,
"gcdout");
1186 TermFree(term1,
"GCDfunction3-a");
1189 MLOCK(ErrorMessageLock);
1190 MesCall(
"GCDfunction3");
1191 MUNLOCK(ErrorMessageLock);
1200 WORD *PutExtraSymbols(
PHEAD WORD *in,WORD startebuf,
int *actionflag)
1202 WORD *termout = AT.WorkPointer;
1211 if ( action > 0 ) *actionflag = 1;
1215 if (
EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 )
goto CalledFrom;
1218 MLOCK(ErrorMessageLock);
1219 MesCall(
"PutExtraSymbols");
1220 MUNLOCK(ErrorMessageLock);
1229 WORD *TakeExtraSymbols(
PHEAD WORD *in,WORD startebuf)
1231 CBUF *C = cbuf+AC.cbufnum;
1232 CBUF *CC = cbuf+AT.ebufnum;
1233 WORD *oldworkpointer = AT.WorkPointer, *termout;
1235 termout = AT.WorkPointer;
1238 if ( ConvertFromPoly(BHEAD in,termout,numxsymbol,CC->numrhs-startebuf+numxsymbol,startebuf-numxsymbol,1) <= 0 ) {
1243 AT.WorkPointer = termout + *termout;
1247 if (
Generator(BHEAD termout,C->numlhs) ) {
1252 AT.WorkPointer = oldworkpointer;
1253 if (
EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 )
goto CalledFrom;
1257 MLOCK(ErrorMessageLock);
1258 MesCall(
"TakeExtraSymbols");
1259 MUNLOCK(ErrorMessageLock);
1268 WORD *MultiplyWithTerm(
PHEAD WORD *in, WORD *term)
1270 WORD *termout, *t, *tt, *tstop, *ttstop;
1271 WORD length, length1, length2;
1272 termout = AT.WorkPointer;
1276 tstop = in + *in; tstop -= ABS(tstop[-1]); t = in + 1;
1277 while ( t < tstop ) *tt++ = *t++;
1278 ttstop = term + *term; ttstop -= ABS(ttstop[-1]); t = term + 1;
1279 while ( t < ttstop ) *tt++ = *t++;
1280 length1 = REDLENG(in[*in-1]); length2 = REDLENG(term[*term-1]);
1281 if ( MulRat(BHEAD (UWORD *)tstop,length1,
1282 (UWORD *)ttstop,length2,(UWORD *)tt,&length) )
goto CalledFrom;
1283 length = INCLENG(length);
1284 tt += ABS(length); tt[-1] = length;
1285 *termout = tt - termout;
1286 Normalize(BHEAD termout);
1290 if (
EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 )
goto CalledFrom;
1294 MLOCK(ErrorMessageLock);
1295 MesCall(
"MultiplyWithTerm");
1296 MUNLOCK(ErrorMessageLock);
1319 WORD *t, *tstop, *tcom, *tout, *tstore, *r, *rstop, *m, *mm, *w, *ww, *wterm;
1320 WORD *tnext, *tt, *tterm, code[2];
1322 int i, j, k, action = 0, sign;
1323 UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMbuffer2, *ap;
1324 WORD GCDlen, GCDlen2, LCMlen, LCMlen2, length, redlength, len1, len2;
1325 tout = tstore = term+1;
1331 tstop = tnext-ABS(tnext[-1]);
1333 while ( t < tstop ) {
1334 if ( *t == INDEX ) {
1335 i = t[1]; NCOPY(tout,t,i);
break;
1339 if ( tout > tstore ) {
1343 rstop = tnext - ABS(tnext[-1]);
1345 while ( r < rstop ) {
1346 if ( *r != INDEX ) { r += r[1];
continue; }
1348 while ( m < tout ) {
1349 for ( i = 2; i < r[1]; i++ ) {
1350 if ( *m == r[i] )
break;
1351 if ( *m > r[i] )
continue;
1353 while ( mm < tout ) { mm[-1] = mm[0]; mm++; }
1354 tout--; tstore[1]--; m--;
1360 if ( r >= rstop || tout <= tstore+2 ) {
1361 tout = tstore;
break;
1364 if ( tout > tstore+2 ) {
1368 tnext = t + *t; t++; w++;
1369 while ( *t != INDEX ) { i = t[1]; NCOPY(w,t,i); }
1370 tt = t + t[1]; t += 2; r = tstore+2; ww = w; *w++ = INDEX; w++;
1371 while ( r < tout && t < tt ) {
1372 if ( *r > *t ) { *w++ = *t++; }
1373 else if ( *r == *t ) { r++; t++; }
1374 else goto CalledFrom;
1376 if ( r < tout )
goto CalledFrom;
1377 while ( t < tt ) *w++ = *t++;
1379 if ( ww[1] == 2 ) w = ww;
1380 while ( t < tnext ) *w++ = *t++;
1391 code[0] = VECTOR; code[1] = DELTA;
1392 for ( k = 0; k < 2; k++ ) {
1395 tstop = tnext-ABS(tnext[-1]);
1397 while ( t < tstop ) {
1398 if ( *t == code[k] ) {
1399 i = t[1]; NCOPY(tout,t,i);
break;
1403 if ( tout > tstore ) {
1407 rstop = tnext - ABS(tnext[-1]);
1409 while ( r < rstop ) {
1410 if ( *r != code[k] ) { r += r[1];
continue; }
1412 while ( m < tout ) {
1413 for ( i = 2; i < r[1]; i += 2 ) {
1414 if ( *m == r[i] && m[1] == r[i+1] )
break;
1415 if ( *m > r[i] || ( *m == r[i] && m[1] > r[i+1] ) )
continue;
1417 while ( mm < tout ) { mm[-2] = mm[0]; mm[-1] = mm[1]; mm += 2; }
1418 tout -= 2; tstore[1] -= 2; m -= 2;
1424 if ( r >= rstop || tout <= tstore+2 ) {
1425 tout = tstore;
break;
1428 if ( tout > tstore+2 ) {
1432 tnext = t + *t; t++; w++;
1433 while ( *t != code[k] ) { i = t[1]; NCOPY(w,t,i); }
1434 tt = t + t[1]; t += 2; r = tstore+2; ww = w; *w++ = code[k]; w++;
1435 while ( r < tout && t < tt ) {
1436 if ( ( *r > *t ) || ( *r == *t && r[1] > t[1] ) )
1437 { *w++ = *t++; *w++ = *t++; }
1438 else if ( *r == *t && r[1] == t[1] ) { r += 2; t += 2; }
1439 else goto CalledFrom;
1441 if ( r < tout )
goto CalledFrom;
1442 while ( t < tt ) *w++ = *t++;
1444 if ( ww[1] == 2 ) w = ww;
1445 while ( t < tnext ) *w++ = *t++;
1459 tstop = tnext-ABS(tnext[-1]);
1462 while ( t < tstop ) {
1463 if ( *t >= FUNCTION ) {
1464 if ( functions[*t-FUNCTION].commute ) {
1465 if ( tcom == 0 ) { tcom = tstore; }
1467 for ( i = 0; i < t[1]; i++ ) {
1468 if ( t[i] != tcom[i] ) {
1469 MLOCK(ErrorMessageLock);
1470 MesPrint(
"GCD or factorization of more than one noncommuting object not allowed");
1471 MUNLOCK(ErrorMessageLock);
1477 i = t[1]; NCOPY(tstore,t,i);
1481 if ( tout > tstore ) {
1484 tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1486 while ( r < tout ) {
1488 while ( tt < tstop ) {
1489 for ( i = 0; i < r[1]; i++ ) {
1490 if ( r[i] != tt[i] )
break;
1492 if ( i == r[1] ) { r += r[1];
goto nextr1; }
1497 m = r; mm = r + r[1];
1498 while ( mm < tout ) *m++ = *mm++;
1502 if ( tout <= tstore )
break;
1506 if ( tout > tstore ) {
1512 while ( r < tout ) {
1513 t = in; ww = in; w = ww+1;
1518 for ( i = 0; i < r[1]; i++ ) {
1519 if ( t[i] != r[i] ) {
1520 j = t[1]; NCOPY(w,t,j);
1526 while ( t < tnext ) *w++ = *t++;
1546 tterm = AT.WorkPointer; tt = tterm+1;
1547 tout[0] = SYMBOL; tout[1] = 2;
1550 tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1551 while ( t < tstop ) {
1552 if ( *t == SYMBOL ) {
1553 MergeSymbolLists(BHEAD tout,t,-1);
1560 if ( tout[1] > 2 ) {
1562 tt[0] = t[0]; tt[1] = t[1];
1563 for ( i = 2; i < t[1]; i += 2 ) {
1564 tt[i] = t[i]; tt[i+1] = -t[i+1];
1579 tout[0] = DOTPRODUCT; tout[1] = 2;
1582 tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1583 while ( t < tstop ) {
1584 if ( *t == DOTPRODUCT ) {
1585 MergeDotproductLists(BHEAD tout,t,-1);
1592 if ( tout[1] > 2 ) {
1594 tt[0] = t[0]; tt[1] = t[1];
1595 for ( i = 2; i < t[1]; i += 3 ) {
1596 tt[i] = t[i]; tt[i+1] = t[i+1]; tt[i+2] = -t[i+2];
1609 AT.WorkPointer = tt;
1610 if ( AN.cmod != 0 ) {
1612 t = in; tnext = t + *t; tstop = tnext - ABS(tnext[-1]);
1614 if ( tnext[-1] < 0 ) x += AC.cmod[0];
1615 if (
GetModInverses(x,(WORD)(AN.cmod[0]),&ix,&ip) )
goto CalledFrom;
1616 *tout++ = x; *tout++ = 1; *tout++ = 3;
1617 *tt++ = ix; *tt++ = 1; *tt++ = 3;
1620 GCDbuffer = NumberMalloc(
"MakeInteger");
1621 GCDbuffer2 = NumberMalloc(
"MakeInteger");
1622 LCMbuffer = NumberMalloc(
"MakeInteger");
1623 LCMbuffer2 = NumberMalloc(
"MakeInteger");
1625 tnext = t + *t; length = tnext[-1];
1626 if ( length < 0 ) { sign = -1; length = -length; }
1628 tstop = tnext - length;
1629 redlength = (length-1)/2;
1630 for ( i = 0; i < redlength; i++ ) {
1631 GCDbuffer[i] = (UWORD)(tstop[i]);
1632 LCMbuffer[i] = (UWORD)(tstop[redlength+i]);
1634 GCDlen = LCMlen = redlength;
1635 while ( GCDbuffer[GCDlen-1] == 0 ) GCDlen--;
1636 while ( LCMbuffer[LCMlen-1] == 0 ) LCMlen--;
1639 tnext = t + *t; length = ABS(tnext[-1]);
1640 tstop = tnext - length; redlength = (length-1)/2;
1641 len1 = len2 = redlength;
1642 den = tstop + redlength;
1643 while ( tstop[len1-1] == 0 ) len1--;
1644 while ( den[len2-1] == 0 ) len2--;
1645 if ( GCDlen == 1 && GCDbuffer[0] == 1 ) {}
1647 GcdLong(BHEAD (UWORD *)tstop,len1,GCDbuffer,GCDlen,GCDbuffer2,&GCDlen2);
1648 ap = GCDbuffer; GCDbuffer = GCDbuffer2; GCDbuffer2 = ap;
1649 a = GCDlen; GCDlen = GCDlen2; GCDlen2 = a;
1651 if ( len2 == 1 && den[0] == 1 ) {}
1653 GcdLong(BHEAD LCMbuffer,LCMlen,(UWORD *)den,len2,LCMbuffer2,&LCMlen2);
1654 DivLong((UWORD *)den,len2,LCMbuffer2,LCMlen2,
1655 GCDbuffer2,&GCDlen2,(UWORD *)AT.WorkPointer,&a);
1656 MulLong(LCMbuffer,LCMlen,GCDbuffer2,GCDlen2,LCMbuffer2,&LCMlen2);
1657 ap = LCMbuffer; LCMbuffer = LCMbuffer2; LCMbuffer2 = ap;
1658 a = LCMlen; LCMlen = LCMlen2; LCMlen2 = a;
1662 if ( GCDlen != 1 || GCDbuffer[0] != 1 || LCMlen != 1 || LCMbuffer[0] != 1 ) {
1663 redlength = GCDlen;
if ( LCMlen > GCDlen ) redlength = LCMlen;
1664 for ( i = 0; i < GCDlen; i++ ) *tout++ = (WORD)(GCDbuffer[i]);
1665 for ( ; i < redlength; i++ ) *tout++ = 0;
1666 for ( i = 0; i < LCMlen; i++ ) *tout++ = (WORD)(LCMbuffer[i]);
1667 for ( ; i < redlength; i++ ) *tout++ = 0;
1668 *tout++ = (2*redlength+1)*sign;
1669 for ( i = 0; i < LCMlen; i++ ) *tt++ = (WORD)(LCMbuffer[i]);
1670 for ( ; i < redlength; i++ ) *tt++ = 0;
1671 for ( i = 0; i < GCDlen; i++ ) *tt++ = (WORD)(GCDbuffer[i]);
1672 for ( ; i < redlength; i++ ) *tt++ = 0;
1673 *tt++ = (2*redlength+1)*sign;
1677 *tout++ = 1; *tout++ = 1; *tout++ = 3*sign;
1678 *tt++ = 1; *tt++ = 1; *tt++ = 3*sign;
1679 if ( sign != 1 ) action++;
1681 NumberFree(LCMbuffer2,
"MakeInteger");
1682 NumberFree(LCMbuffer ,
"MakeInteger");
1683 NumberFree(GCDbuffer2,
"MakeInteger");
1684 NumberFree(GCDbuffer ,
"MakeInteger");
1691 *tterm = tt - tterm;
1692 AT.WorkPointer = tt;
1693 inp = MultiplyWithTerm(BHEAD in,tterm);
1694 AT.WorkPointer = tterm;
1700 *term = tout - term;
1701 AT.WorkPointer = tterm;
1704 MLOCK(ErrorMessageLock);
1705 MesCall(
"TakeContent");
1706 MUNLOCK(ErrorMessageLock);
1722 int MergeSymbolLists(
PHEAD WORD *old, WORD *extra,
int par)
1725 WORD *
new = TermMalloc(
"MergeSymbolLists");
1726 WORD *t1, *t2, *fill;
1729 i1 = old[1] - 2; i2 = extra[1] - 2;
1730 t1 = old + 2; t2 = extra + 2;
1733 while ( i1 > 0 && i2 > 0 ) {
1735 if ( t2[1] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1738 else if ( *t1 < *t2 ) {
1739 if ( t1[1] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1742 else if ( t1[1] < t2[1] ) {
1743 *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1746 *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1751 while ( i1 > 0 && i2 > 0 ) {
1753 if ( t2[1] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1756 else if ( *t1 < *t2 ) {
1757 if ( t1[1] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1760 else if ( t1[1] > t2[1] ) {
1761 *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1764 *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1769 while ( i1 > 0 && i2 > 0 ) {
1773 else if ( *t1 < *t2 ) {
1776 else if ( ( t1[1] > 0 ) && ( t2[1] < 0 ) ) { t1 += 2; t2 += 2; }
1777 else if ( ( t1[1] < 0 ) && ( t2[1] > 0 ) ) { t1 += 2; t2 += 2; }
1778 else if ( t1[1] > 0 ) {
1779 if ( t1[1] < t2[1] ) {
1780 *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1783 *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1787 if ( t2[1] < t1[1] ) {
1788 *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1791 *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1797 i1 =
new[1] = fill -
new;
1798 t2 =
new; t1 = old; NCOPY(t1,t2,i1);
1799 TermFree(
new,
"MergeSymbolLists");
1815 int MergeDotproductLists(
PHEAD WORD *old, WORD *extra,
int par)
1818 WORD *
new = TermMalloc(
"MergeDotproductLists");
1819 WORD *t1, *t2, *fill;
1822 i1 = old[1] - 2; i2 = extra[1] - 2;
1823 t1 = old + 2; t2 = extra + 2;
1826 while ( i1 > 0 && i2 > 0 ) {
1827 if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1828 if ( t2[2] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; }
1831 else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1832 if ( t1[2] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; }
1835 else if ( t1[2] < t2[2] ) {
1836 *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1839 *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1844 while ( i1 > 0 && i2 > 0 ) {
1845 if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1846 if ( t2[2] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; }
1849 else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1850 if ( t1[2] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; }
1853 else if ( t1[2] > t2[2] ) {
1854 *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1857 *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1862 while ( i1 > 0 && i2 > 0 ) {
1863 if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1866 else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1869 else if ( ( t1[2] > 0 ) && ( t2[2] < 0 ) ) { t1 += 3; t2 += 3; }
1870 else if ( ( t1[2] < 0 ) && ( t2[2] > 0 ) ) { t1 += 3; t2 += 3; }
1871 else if ( t1[2] > 0 ) {
1872 if ( t1[2] < t2[2] ) {
1873 *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1876 *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1880 if ( t2[2] < t1[2] ) {
1881 *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1884 *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1890 i1 =
new[1] = fill -
new;
1891 t2 =
new; t1 = old; NCOPY(t1,t2,i1);
1892 TermFree(
new,
"MergeDotproductLists");
1908 WORD *CreateExpression(
PHEAD WORD nexp)
1911 CBUF *C = cbuf+AC.cbufnum;
1912 POSITION startposition, oldposition;
1914 WORD *term, *oldipointer = AR.CompressPointer;
1916 switch ( Expressions[nexp].status ) {
1917 case HIDDENLEXPRESSION:
1918 case HIDDENGEXPRESSION:
1919 case DROPHLEXPRESSION:
1920 case DROPHGEXPRESSION:
1921 case UNHIDELEXPRESSION:
1922 case UNHIDEGEXPRESSION:
1923 AR.GetOneFile = 2; fi = AR.hidefile;
1926 AR.GetOneFile = 0; fi = AR.infile;
1929 SeekScratch(fi,&oldposition);
1930 startposition = AS.OldOnFile[nexp];
1931 term = AT.WorkPointer;
1932 if ( GetOneTerm(BHEAD term,fi,&startposition,0) <= 0 )
goto CalledFrom;
1934 AR.CompressPointer = oldipointer;
1935 while ( GetOneTerm(BHEAD term,fi,&startposition,0) > 0 ) {
1936 AT.WorkPointer = term + *term;
1937 if (
Generator(BHEAD term,C->numlhs) ) {
1941 AR.CompressPointer = oldipointer;
1943 AT.WorkPointer = term;
1944 if (
EndSort(BHEAD (WORD *)((VOID *)(&term)),2) < 0 )
goto CalledFrom;
1945 SetScratch(fi,&oldposition);
1948 MLOCK(ErrorMessageLock);
1949 MesCall(
"CreateExpression");
1950 MUNLOCK(ErrorMessageLock);
1964 int GCDterms(
PHEAD WORD *term1, WORD *term2, WORD *termout)
1967 WORD *t1, *t1stop, *t1next, *t2, *t2stop, *t2next, *tout, *tt1, *tt2;
1968 int count1, count2, i, ii, x1, sign;
1969 WORD length1, length2;
1970 t1 = term1 + *term1; t1stop = t1 - ABS(t1[-1]); t1 = term1+1;
1971 t2 = term2 + *term2; t2stop = t2 - ABS(t2[-1]); t2 = term2+1;
1973 while ( t1 < t1stop ) {
1974 t1next = t1 + t1[1];
1976 if ( *t1 == SYMBOL ) {
1977 while ( t2 < t2stop && *t2 != SYMBOL ) t2 += t2[1];
1978 if ( *t2 == SYMBOL ) {
1980 tt1 = t1+2; tt2 = t2+2; count1 = 0;
1981 while ( tt1 < t1next && tt2 < t2next ) {
1982 if ( *tt1 < *tt2 ) tt1 += 2;
1983 else if ( *tt1 > *tt2 ) tt2 += 2;
1984 else if ( ( tt1[1] > 0 && tt2[1] < 0 ) ||
1985 ( tt2[1] > 0 && tt1[1] < 0 ) ) {
1990 if ( tt1[1] < 0 ) {
if ( tt2[1] > x1 ) x1 = tt2[1]; }
1991 else {
if ( tt2[1] < x1 ) x1 = tt2[1]; }
1992 tout[count1+2] = *tt1;
1993 tout[count1+3] = x1;
1999 *tout = SYMBOL; tout[1] = count1+2; tout += tout[1];
2003 else if ( *t1 == DOTPRODUCT ) {
2004 while ( t2 < t2stop && *t2 != DOTPRODUCT ) t2 += t2[1];
2005 if ( *t2 == DOTPRODUCT ) {
2007 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2008 while ( tt1 < t1next && tt2 < t2next ) {
2009 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 3;
2010 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 3;
2011 else if ( ( tt1[2] > 0 && tt2[2] < 0 ) ||
2012 ( tt2[2] > 0 && tt1[2] < 0 ) ) {
2017 if ( tt1[2] < 0 ) {
if ( tt2[2] > x1 ) x1 = tt2[2]; }
2018 else {
if ( tt2[2] < x1 ) x1 = tt2[2]; }
2019 tout[count1+2] = *tt1;
2020 tout[count1+3] = tt1[1];
2021 tout[count1+4] = x1;
2027 *tout = DOTPRODUCT; tout[1] = count1+2; tout += tout[1];
2031 else if ( *t1 == VECTOR ) {
2032 while ( t2 < t2stop && *t2 != VECTOR ) t2 += t2[1];
2033 if ( *t2 == VECTOR ) {
2035 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2036 while ( tt1 < t1next && tt2 < t2next ) {
2037 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2038 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2040 tout[count1+2] = *tt1;
2041 tout[count1+3] = tt1[1];
2047 *tout = VECTOR; tout[1] = count1+2; tout += tout[1];
2051 else if ( *t1 == INDEX ) {
2052 while ( t2 < t2stop && *t2 != INDEX ) t2 += t2[1];
2053 if ( *t2 == INDEX ) {
2055 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2056 while ( tt1 < t1next && tt2 < t2next ) {
2057 if ( *tt1 < *tt2 ) tt1 += 1;
2058 else if ( *tt1 > *tt2 ) tt2 += 1;
2060 tout[count1+2] = *tt1;
2066 *tout = INDEX; tout[1] = count1+2; tout += tout[1];
2070 else if ( *t1 == DELTA ) {
2071 while ( t2 < t2stop && *t2 != DELTA ) t2 += t2[1];
2072 if ( *t2 == DELTA ) {
2074 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2075 while ( tt1 < t1next && tt2 < t2next ) {
2076 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2077 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2079 tout[count1+2] = *tt1;
2080 tout[count1+3] = tt1[1];
2086 *tout = DELTA; tout[1] = count1+2; tout += tout[1];
2090 else if ( *t1 >= FUNCTION ) {
2096 while ( t1next < t1stop && *t1 == *t1next && t1[1] == t1next[1] ) {
2097 for ( i = 2; i < t1[1]; i++ ) {
2098 if ( t1[i] != t1next[i] )
break;
2100 if ( i < t1[1] )
break;
2102 t1next += t1next[1];
2105 while ( t2 < t2stop ) {
2106 if ( *t2 == *t1 && t2[1] == t1[1] ) {
2107 for ( i = 2; i < t1[1]; i++ ) {
2108 if ( t2[i] != t1[i] )
break;
2110 if ( i >= t1[1] ) count2++;
2114 if ( count1 < count2 ) count2 = count1;
2117 while ( count2 > 0 ) { tout += tout[1]; count2--; }
2131 length1 = term1[*term1-1]; ii = i = ABS(length1); t1 = t1stop;
2132 if ( t1 != tout ) { NCOPY(tout,t1,i); tout -= ii; }
2133 length2 = term2[*term2-1];
2134 if ( length1 < 0 && length2 < 0 ) sign = -1;
2135 if ( AccumGCD(BHEAD (UWORD *)tout,&length1,(UWORD *)t2stop,length2) ) {
2136 MLOCK(ErrorMessageLock);
2137 MesCall(
"GCDterms");
2138 MUNLOCK(ErrorMessageLock);
2141 if ( sign < 0 && length1 > 0 ) length1 = -length1;
2142 tout += ABS(length1); tout[-1] = length1;
2143 *termout = tout - termout;
2161 WORD divrem[3] = { DIVFUNCTION, REMFUNCTION, INVERSEFUNCTION };
2163 int DIVfunction(
PHEAD WORD *term,WORD level,
int par)
2166 WORD *t, *tt, *r, *arg1 = 0, *arg2 = 0, *arg3 = 0, *termout;
2167 WORD *tstop, *tend, *r3, *rr, *rstop, tlength, rlength, newlength;
2168 WORD *proper1, *proper2, *proper3 = 0;
2169 int numargs = 0, type1, type2, actionflag1, actionflag2;
2170 WORD startebuf = cbuf[AT.ebufnum].numrhs;
2171 if ( par < 0 || par > 2 ) {
2172 MLOCK(ErrorMessageLock);
2173 MesPrint(
"Internal error. Illegal parameter %d in DIVfunction.",par);
2174 MUNLOCK(ErrorMessageLock);
2180 tend = term + *term; tstop = tend - ABS(tend[-1]);
2182 while ( t < tstop ) {
2183 if ( *t != divrem[par] ) { t += t[1];
continue; }
2185 tt = t + t[1]; numargs = 0;
2187 if ( numargs == 0 ) { arg1 = r; }
2188 if ( numargs == 1 ) { arg2 = r; }
2192 if ( numargs == 2 )
break;
2196 MLOCK(ErrorMessageLock);
2197 MesPrint(
"Internal error. Indicated div_ or rem_ function not encountered.");
2198 MUNLOCK(ErrorMessageLock);
2204 if ( *arg1 == -SNUMBER && arg1[1] == 0 ) {
2205 if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2207 MLOCK(ErrorMessageLock);
2208 MesPrint(
"0/0 in either div_ or rem_ function.");
2209 MUNLOCK(ErrorMessageLock);
2214 if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2216 MLOCK(ErrorMessageLock);
2217 MesPrint(
"Division by zero in either div_ or rem_ function.");
2218 MUNLOCK(ErrorMessageLock);
2221 if ( ( arg1 = ConvertArgument(BHEAD arg1, &type1) ) == 0 )
goto CalledFrom;
2222 if ( ( arg2 = ConvertArgument(BHEAD arg2, &type2) ) == 0 )
goto CalledFrom;
2225 M_free(arg2,
"DIVfunction");
2226 M_free(arg1,
"DIVfunction");
2229 M_free(arg2,
"DIVfunction");
2230 M_free(arg1,
"DIVfunction");
2234 M_free(arg2,
"DIVfunction");
2235 M_free(arg1,
"DIVfunction");
2238 if ( ( proper1 = PutExtraSymbols(BHEAD arg1,startebuf,&actionflag1) ) == 0 )
goto CalledFrom;
2239 if ( ( proper2 = PutExtraSymbols(BHEAD arg2,startebuf,&actionflag2) ) == 0 )
goto CalledFrom;
2248 M_free(arg2,
"DIVfunction");
2257 M_free(arg1,
"DIVfunction");
2258 if ( par == 0 ) proper3 = poly_div(BHEAD proper1, proper2);
2259 else if ( par == 1 ) proper3 = poly_rem(BHEAD proper1, proper2);
2260 else if ( par == 2 ) proper3 = poly_inverse(BHEAD proper1, proper2);
2261 if ( proper3 == 0 )
goto CalledFrom;
2262 if ( actionflag1 || actionflag2 ) {
2263 if ( ( arg3 = TakeExtraSymbols(BHEAD proper3,startebuf) ) == 0 )
goto CalledFrom;
2264 M_free(proper3,
"DIVfunction");
2269 M_free(proper2,
"DIVfunction");
2270 M_free(proper1,
"DIVfunction");
2271 cbuf[AT.ebufnum].numrhs = startebuf;
2273 termout = AT.WorkPointer;
2275 tlength = REDLENG(tlength);
2278 tt = term + 1; rr = termout + 1;
2279 while ( tt < t ) *rr++ = *tt++;
2282 rstop = r3 - ABS(r3[-1]);
2283 while ( r < rstop ) *rr++ = *r++;
2285 while ( tt < tstop ) *rr++ = *tt++;
2287 rlength = REDLENG(rlength);
2288 if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)rstop,rlength,
2289 (UWORD *)rr,&newlength) < 0 )
goto CalledFrom;
2290 rlength = INCLENG(newlength);
2293 *termout = rr - termout;
2294 AT.WorkPointer = rr;
2295 if (
Generator(BHEAD termout,level) )
goto CalledFrom;
2297 AT.WorkPointer = termout;
2299 M_free(arg3,
"DIVfunction");
2302 MLOCK(ErrorMessageLock);
2303 MesCall(
"DIVfunction");
2304 MUNLOCK(ErrorMessageLock);
2315 WORD *ConvertArgument(
PHEAD WORD *arg,
int *type)
2317 WORD *output, *t, *r;
2320 output = (WORD *)Malloc1((*arg)*
sizeof(WORD),
"ConvertArgument");
2321 i = *arg - ARGHEAD; t = arg + ARGHEAD; r = output;
2326 if ( *arg == -EXPRESSION ) {
2328 return(CreateExpression(BHEAD arg[1]));
2330 if ( *arg == -DOLLAREXPRESSION ) {
2333 d = DolToTerms(BHEAD arg[1]);
2339 output = (WORD *)Malloc1((d->size+1)*
sizeof(WORD),
"Copy of dollar content");
2340 WCOPY(output,d->where,d->size+1);
2341 if ( d->factors ) { M_free(d->factors,
"Dollar factors"); d->factors = 0; }
2342 M_free(d,
"Copy of dollar variable");
2350 output = (WORD *)Malloc1(size*
sizeof(WORD),
"ConvertArgument");
2353 output[0] = 8; output[1] = SYMBOL; output[2] = 4; output[3] = arg[1];
2354 output[4] = 1; output[5] = 1; output[6] = 1; output[7] = 3;
2360 output[0] = 7; output[1] = INDEX; output[2] = 3; output[3] = arg[1];
2361 output[4] = 1; output[5] = 1;
2362 if ( *arg == -MINVECTOR ) output[6] = -3;
2369 output[1] = -arg[1]; output[2] = 1; output[3] = -3;
2372 output[1] = arg[1]; output[2] = 1; output[3] = 3;
2377 output[0] = FUNHEAD+4;
2379 output[2] = FUNHEAD;
2380 for ( i = 3; i <= FUNHEAD; i++ ) output[i] = 0;
2381 output[FUNHEAD+1] = 1;
2382 output[FUNHEAD+2] = 1;
2383 output[FUNHEAD+3] = 3;
2384 output[FUNHEAD+4] = 0;
WORD * TakeContent(PHEAD WORD *in, WORD *term)
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
WORD InsertTerm(PHEAD WORD *, WORD, WORD, WORD *, WORD *, WORD)
int GetModInverses(WORD, WORD, WORD *, WORD *)
WORD StoreTerm(PHEAD WORD *)
WORD Generator(PHEAD WORD *, WORD)
LONG EndSort(PHEAD WORD *, int)