48 WORD CompareFunctions(WORD *fleft,WORD *fright)
51 if ( AC.properorderflag ) {
52 if ( ( *fleft >= (FUNCTION+WILDOFFSET)
53 && functions[*fleft-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
54 || ( *fleft >= FUNCTION && *fleft < (FUNCTION + WILDOFFSET)
55 && functions[*fleft-FUNCTION].spec >= TENSORFUNCTION ) ) {}
57 WORD *s1, *s2, *ss1, *ss2;
58 s1 = fleft+FUNHEAD; s2 = fright+FUNHEAD;
59 ss1 = fleft + fleft[1]; ss2 = fright + fright[1];
60 while ( s1 < ss1 && s2 < ss2 ) {
62 if ( k > 0 )
return(1);
63 if ( k < 0 )
return(0);
67 if ( s1 < ss1 )
return(1);
70 k = fleft[1] - FUNHEAD;
71 kk = fright[1] - FUNHEAD;
74 while ( k > 0 && kk > 0 ) {
75 if ( *fleft < *fright )
return(0);
76 else if ( *fleft++ > *fright++ )
return(1);
79 if ( k > 0 )
return(1);
83 k = fleft[1] - FUNHEAD;
84 kk = fright[1] - FUNHEAD;
87 while ( k > 0 && kk > 0 ) {
88 if ( *fleft < *fright )
return(0);
89 else if ( *fleft++ > *fright++ )
return(1);
92 if ( k > 0 )
return(1);
112 WORD Commute(WORD *fleft, WORD *fright)
115 if ( *fleft == DOLLAREXPRESSION || *fright == DOLLAREXPRESSION )
return(0);
116 fun1 = ABS(*fleft); fun2 = ABS(*fright);
117 if ( *fleft >= GAMMA && *fleft <= GAMMASEVEN
118 && *fright >= GAMMA && *fright <= GAMMASEVEN ) {
119 if ( fleft[FUNHEAD] < AM.OffsetIndex && fleft[FUNHEAD] > fright[FUNHEAD] )
123 if ( fun1 >= WILDOFFSET ) fun1 -= WILDOFFSET;
124 if ( fun2 >= WILDOFFSET ) fun2 -= WILDOFFSET;
125 if ( ( ( functions[fun1-FUNCTION].flags & COULDCOMMUTE ) == 0 )
126 || ( ( functions[fun2-FUNCTION].flags & COULDCOMMUTE ) == 0 ) )
return(0);
131 if ( AC.CommuteInSet == 0 )
return(0);
146 if ( fun1 >= fun2 ) {
147 WORD *group = AC.CommuteInSet, *g1, *g2, *g3;
148 while ( *group > 0 ) {
152 if ( *g1 == fun1 || ( fun1 <= GAMMASEVEN && fun1 >= GAMMA
153 && *g1 <= GAMMASEVEN && *g1 >= GAMMA ) ) {
156 if ( g1 != g2 && ( *g2 == fun2 ||
157 ( fun2 <= GAMMASEVEN && fun2 >= GAMMA
158 && *g2 <= GAMMASEVEN && *g2 >= GAMMA ) ) ) {
159 if ( fun1 != fun2 )
return(1);
160 if ( *fleft < 0 )
return(0);
161 if ( *fright < 0 )
return(1);
162 return(CompareFunctions(fleft,fright));
187 WORD Normalize(
PHEAD WORD *term)
193 WORD *t, *m, *r, i, j, k, l, nsym, *ss, *tt, *u;
194 WORD shortnum, stype;
195 WORD *stop, *to = 0, *from = 0;
201 WORD psym[7*NORMSIZE],*ppsym;
202 WORD pvec[NORMSIZE],*ppvec,nvec;
203 WORD pdot[3*NORMSIZE],*ppdot,ndot;
204 WORD pdel[2*NORMSIZE],*ppdel,ndel;
205 WORD pind[NORMSIZE],nind;
206 WORD *peps[NORMSIZE/3],neps;
207 WORD *pden[NORMSIZE/3],nden;
208 WORD *pcom[NORMSIZE],ncom;
209 WORD *pnco[NORMSIZE],nnco;
210 WORD *pcon[2*NORMSIZE],ncon;
212 WORD *n_llnum, *lnum, nnum;
213 WORD *termout, oldtoprhs = 0, subtype;
214 WORD ReplaceType, ReplaceVeto = 0, didcontr, regval = 0;
217 CBUF *C = cbuf+AT.ebufnum;
218 WORD *ANsc = 0, *ANsm = 0, *ANsr = 0, PolyFunMode;
219 LONG oldcpointer = 0;
220 n_coef = TermMalloc(
"NormCoef");
221 n_llnum = TermMalloc(
"n_llnum");
237 if ( !*t ) { TermFree(n_coef,
"NormCoef"); TermFree(n_llnum,
"n_llnum");
return(regval); }
245 termout = AT.WorkPointer;
246 AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
247 fillsetexp = termout+1;
248 AN.PolyNormFlag = 0; PolyFunMode = AN.PolyFunTodo;
256 nsym = nvec = ndot = ndel = neps = nden =
257 nind = ncom = nnco = ncon = 0;
276 if ( *t <= DENOMINATORSYMBOL && *t >= COEFFSYMBOL ) {
282 if ( AN.cTerm ) m = AN.cTerm;
285 ncoef = REDLENG(ncoef);
286 if ( *t == COEFFSYMBOL ) {
288 nnum = REDLENG(m[-1]);
292 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,nnum,
293 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
299 if ( DivRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,nnum,
300 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
308 if ( *t == NUMERATORSYMBOL ) { m -= nnum + 1; }
310 while ( *m == 0 && nnum > 1 ) { m--; nnum--; }
312 if ( i < 0 && *t == NUMERATORSYMBOL ) nnum = -nnum;
316 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)m,nnum) )
323 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)m,nnum) )
329 ncoef = INCLENG(ncoef);
333 else if ( *t == DIMENSIONSYMBOL ) {
334 if ( AN.cTerm ) m = AN.cTerm;
336 k = DimensionTerm(m);
337 if ( k == 0 )
goto NormZero;
338 if ( k == MAXPOSITIVE ) {
339 MLOCK(ErrorMessageLock);
340 MesPrint(
"Dimension_ is undefined in term %t");
341 MUNLOCK(ErrorMessageLock);
344 if ( k == -MAXPOSITIVE ) {
345 MLOCK(ErrorMessageLock);
346 MesPrint(
"Dimension_ out of range in term %t");
347 MUNLOCK(ErrorMessageLock);
350 if ( k > 0 ) { *((UWORD *)lnum) = k; nnum = 1; }
351 else { *((UWORD *)lnum) = -k; nnum = -1; }
352 ncoef = REDLENG(ncoef);
353 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
goto FromNorm;
354 ncoef = INCLENG(ncoef);
358 if ( ( *t >= MAXPOWER && *t < 2*MAXPOWER )
359 || ( *t < -MAXPOWER && *t > -2*MAXPOWER ) ) {
366 if ( t[1] & 1 ) ncoef = -ncoef;
368 else if ( *t == MAXPOWER ) {
369 if ( t[1] > 0 )
goto NormZero;
377 if ( t[1] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[1]))) )
379 ncoef = REDLENG(ncoef);
381 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
384 else if ( t[1] > 0 ) {
385 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
388 ncoef = INCLENG(ncoef);
395 if ( ( *t <= NumSymbols && *t > -MAXPOWER )
396 && ( symbols[*t].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
397 if ( t[1] <= 2*MAXPOWER && t[1] >= -2*MAXPOWER ) {
398 t[1] %= symbols[*t].maxpower;
399 if ( t[1] < 0 ) t[1] += symbols[*t].maxpower;
400 if ( ( symbols[*t].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
401 if ( ( ( symbols[*t].maxpower & 1 ) == 0 ) &&
402 ( t[1] >= symbols[*t].maxpower/2 ) ) {
403 t[1] -= symbols[*t].maxpower/2; ncoef = -ncoef;
406 if ( t[1] == 0 ) { t += 2;
goto NextSymbol; }
415 if ( *t > 2*MAXPOWER || *t < -2*MAXPOWER
416 || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
417 MLOCK(ErrorMessageLock);
418 MesPrint(
"Illegal wildcard power combination.");
419 MUNLOCK(ErrorMessageLock);
423 if ( ( t[-1] <= NumSymbols && t[-1] > -MAXPOWER )
424 && ( symbols[t[-1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
425 *m %= symbols[t[-1]].maxpower;
426 if ( *m < 0 ) *m += symbols[t[-1]].maxpower;
427 if ( ( symbols[t[-1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
428 if ( ( ( symbols[t[-1]].maxpower & 1 ) == 0 ) &&
429 ( *m >= symbols[t[-1]].maxpower/2 ) ) {
430 *m -= symbols[t[-1]].maxpower/2; ncoef = -ncoef;
434 if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
435 MLOCK(ErrorMessageLock);
436 MesPrint(
"Power overflow during normalization");
437 MUNLOCK(ErrorMessageLock);
443 { *m = m[2]; m++; *m = m[2]; m++; i++; }
450 }
while ( *t < *m && --i > 0 );
453 { m--; m[2] = *m; m--; m[2] = *m; i++; }
465 if ( t[1] == FUNNYVEC ) {
469 else if ( t[1] < 0 ) {
470 if ( *t == NOINDEX && t[1] == NOINDEX ) t += 2;
472 *ppdot++ = *t++; *ppdot++ = *t++; *ppdot++ = 1; ndot++;
475 else { *ppvec++ = *t++; *ppvec++ = *t++; nvec += 2; }
481 if ( t[2] == 0 ) t += 3;
482 else if ( ndot > 0 && t[0] == ppdot[-3]
483 && t[1] == ppdot[-2] ) {
486 if ( ppdot[-1] == 0 ) { ppdot -= 3; ndot--; }
489 *ppdot++ = *t++; *ppdot++ = *t++;
490 *ppdot++ = *t++; ndot++;
497 if ( WildFill(BHEAD termout,term,AT.dummysubexp) < 0 )
goto FromNorm;
499 t = termout; m = term;
502 case DOLLAREXPRESSION :
509 if ( AR.Eside != LHSIDE ) {
512 int nummodopt, ptype = -1;
513 if ( AS.MultiThreaded ) {
514 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
515 if ( t[2] == ModOptdollars[nummodopt].number )
break;
517 if ( nummodopt < NumModOptdollars ) {
518 ptype = ModOptdollars[nummodopt].type;
519 if ( ptype == MODLOCAL ) {
520 d = ModOptdollars[nummodopt].dstruct+AT.identity;
523 LOCK(d->pthreadslockread);
528 if ( d->type == DOLZERO ) {
530 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
532 if ( t[3] == 0 )
goto NormZZ;
533 if ( t[3] < 0 )
goto NormInf;
536 else if ( d->type == DOLNUMBER ) {
539 nnum = d->where[nnum-1];
540 if ( nnum < 0 ) { ncoef = -ncoef; nnum = -nnum; }
542 for ( i = 1; i <= nnum; i++ ) lnum[i-1] = d->where[i];
544 if ( nnum == 0 || ( nnum == 1 && lnum[0] == 0 ) ) {
546 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
548 if ( t[3] < 0 )
goto NormInf;
549 else if ( t[3] == 0 )
goto NormZZ;
552 if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) )
goto FromNorm;
553 ncoef = REDLENG(ncoef);
555 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) {
557 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
562 else if ( t[3] > 0 ) {
563 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) {
565 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
570 ncoef = INCLENG(ncoef);
572 else if ( d->type == DOLINDEX ) {
573 if ( d->index == 0 ) {
575 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
579 if ( d->index != NOINDEX ) pind[nind++] = d->index;
581 else if ( d->type == DOLTERMS ) {
582 if ( t[3] >= MAXPOWER || t[3] <= -MAXPOWER ) {
583 if ( d->where[0] == 0 )
goto NormZero;
584 if ( d->where[d->where[0]] != 0 ) {
586 MLOCK(ErrorMessageLock);
587 MesPrint(
"!!!Illegal $ expansion with wildcard power!!!");
588 MUNLOCK(ErrorMessageLock);
595 { WORD *td, *tdstop, dj;
597 tdstop = d->where+d->where[0];
598 if ( tdstop[-1] != 3 || tdstop[-2] != 1
599 || tdstop[-3] != 1 )
goto IllDollarExp;
601 if ( td >= tdstop )
goto IllDollarExp;
602 while ( td < tdstop ) {
603 if ( *td == SYMBOL ) {
604 for ( dj = 2; dj < td[1]; dj += 2 ) {
605 if ( td[dj+1] == 1 ) {
610 else if ( td[dj+1] == -1 ) {
615 else goto IllDollarExp;
618 else if ( *td == DOTPRODUCT ) {
619 for ( dj = 2; dj < td[1]; dj += 3 ) {
620 if ( td[dj+2] == 1 ) {
626 else if ( td[dj+2] == -1 ) {
632 else goto IllDollarExp;
635 else goto IllDollarExp;
642 t[0] = SUBEXPRESSION;
646 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
653 if ( *t == DOLLAREXPRESSION ) {
655 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
659 if ( AS.MultiThreaded ) {
660 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
661 if ( t[2] == ModOptdollars[nummodopt].number )
break;
663 if ( nummodopt < NumModOptdollars ) {
664 ptype = ModOptdollars[nummodopt].type;
665 if ( ptype == MODLOCAL ) {
666 d = ModOptdollars[nummodopt].dstruct+AT.identity;
669 LOCK(d->pthreadslockread);
674 if ( d->type == DOLTERMS ) {
675 t[0] = SUBEXPRESSION;
682 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
688 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
690 MLOCK(ErrorMessageLock);
691 MesPrint(
"!!!This $ variation has not been implemented yet!!!");
692 MUNLOCK(ErrorMessageLock);
696 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
712 if ( *t == SUMMEDIND ) {
713 if ( t[1] < -NMIN4SHIFT ) {
714 k = -t[1]-NMIN4SHIFT;
715 k = ExtraSymbol(k,1,nsym,ppsym,&ncoef);
719 else if ( t[1] == 0 )
goto NormZero;
729 ncoef = REDLENG(ncoef);
730 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
732 ncoef = INCLENG(ncoef);
736 else if ( *t == NOINDEX && t[1] == NOINDEX ) t += 2;
737 else if ( *t == EMPTYINDEX && t[1] == EMPTYINDEX ) {
738 *ppdel++ = *t++; *ppdel++ = *t++; ndel += 2;
742 *ppdot++ = *t++; *ppdot++ = *t++; *ppdot++ = 1; ndot++;
745 *ppvec++ = *t++; *ppvec++ = *t++; nvec += 2;
750 *ppvec++ = t[1]; *ppvec++ = *t; t+=2; nvec += 2;
752 else { *ppdel++ = *t++; *ppdel++ = *t++; ndel += 2; }
760 if ( t[FUNHEAD] == -SNUMBER && t[1] == FUNHEAD+2
761 && t[FUNHEAD+1] >= 0 ) {
762 if ( Factorial(BHEAD t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
765 ncoef = REDLENG(ncoef);
766 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
goto FromNorm;
767 ncoef = INCLENG(ncoef);
769 else pcom[ncom++] = t;
771 case BERNOULLIFUNCTION :
775 if ( ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] >= 0 )
776 && ( t[1] == FUNHEAD+2 || ( t[1] == FUNHEAD+4 &&
777 t[FUNHEAD+2] == -SNUMBER && ABS(t[FUNHEAD+3]) == 1 ) ) ) {
779 if ( Bernoulli(t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
781 if ( nnum == 0 )
goto NormZero;
782 inum = nnum;
if ( inum < 0 ) inum = -inum;
785 while ( lnum[mnum-1] == 0 ) mnum--;
786 if ( nnum < 0 ) mnum = -mnum;
787 ncoef = REDLENG(ncoef);
788 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,mnum) )
goto FromNorm;
790 while ( lnum[inum+mnum-1] == 0 ) mnum--;
791 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(lnum+inum),mnum) )
goto FromNorm;
792 ncoef = INCLENG(ncoef);
793 if ( t[1] == FUNHEAD+4 && t[FUNHEAD+1] == 1
794 && t[FUNHEAD+3] == -1 ) ncoef = -ncoef;
796 else pcom[ncom++] = t;
808 if ( k == 0 )
goto NormZero;
809 *((UWORD *)lnum) = k;
817 if ( *t == -EXPRESSION ) {
818 k = AS.OldNumFactors[t[1]];
820 else if ( *t == -DOLLAREXPRESSION ) {
821 k = Dollars[t[1]].nfactors;
827 if ( k == 0 )
goto NormZero;
828 *((UWORD *)lnum) = k;
835 if ( t[FUNHEAD] < 0 ) {
836 if ( t[FUNHEAD] <= -FUNCTION && t[1] == FUNHEAD+1 )
break;
837 if ( t[FUNHEAD] > -FUNCTION && t[1] == FUNHEAD+2 ) {
838 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] == 0 )
goto NormZero;
844 if ( t[FUNHEAD] > 0 && t[FUNHEAD] == t[1]-FUNHEAD ) {
846 t += FUNHEAD+ARGHEAD;
851 if ( k == 0 )
goto NormZero;
852 *((UWORD *)lnum) = k;
856 else pcom[ncom++] = t;
860 k = CountFun(AN.cTerm,t);
863 k = CountFun(term,t);
865 if ( k == 0 )
goto NormZero;
866 if ( k > 0 ) { *((UWORD *)lnum) = k; nnum = 1; }
867 else { *((UWORD *)lnum) = -k; nnum = -1; }
871 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+2] == -SNUMBER
872 && t[1] == FUNHEAD+4 && t[FUNHEAD+3] > 1 ) {
874 if ( t[FUNHEAD+1] == 0 )
goto NormZero;
875 if ( t[FUNHEAD+1] < 0 ) { t[FUNHEAD+1] = -t[FUNHEAD+1]; sgn = -1; }
877 if ( MakeRational(t[FUNHEAD+1],t[FUNHEAD+3],x1,x1+1) ) {
878 static int warnflag = 1;
880 MesPrint(
"%w Warning: fraction could not be reconstructed in MakeRational_");
883 x1[0] = t[FUNHEAD+1]; x1[1] = 1;
885 if ( sgn < 0 ) { t[FUNHEAD+1] = -t[FUNHEAD+1]; x1[0] = -x1[0]; }
886 if ( x1[0] < 0 ) { sgn = -1; x1[0] = -x1[0]; }
888 ncoef = REDLENG(ncoef);
889 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)x1,sgn,
890 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
891 ncoef = INCLENG(ncoef);
894 WORD narg = 0, *tt, *ttstop, *arg1 = 0, *arg2 = 0;
897 ttstop = t + t[1]; tt = t+FUNHEAD;
898 while ( tt < ttstop ) {
900 if ( narg == 1 ) arg1 = tt;
904 if ( narg != 2 )
goto defaultcase;
905 if ( *arg2 == -SNUMBER && arg2[1] <= 1 )
goto defaultcase;
906 else if ( *arg2 > 0 && ttstop[-1] < 0 )
goto defaultcase;
907 x1 = NumberMalloc(
"Norm-MakeRational");
908 if ( *arg1 == -SNUMBER ) {
909 if ( arg1[1] == 0 ) {
910 NumberFree(x1,
"Norm-MakeRational");
913 if ( arg1[1] < 0 ) { x1[0] = -arg1[1]; nx1 = -1; }
914 else { x1[0] = arg1[1]; nx1 = 1; }
916 else if ( *arg1 > 0 ) {
918 nx1 = (ABS(arg2[-1])-1)/2;
919 tc = arg1+ARGHEAD+1+nx1;
921 NumberFree(x1,
"Norm-MakeRational");
924 for ( i = 1; i < nx1; i++ )
if ( tc[i] != 0 ) {
925 NumberFree(x1,
"Norm-MakeRational");
929 for ( i = 0; i < nx1; i++ ) x1[i] = tc[i];
930 if ( arg2[-1] < 0 ) nx1 = -nx1;
933 NumberFree(x1,
"Norm-MakeRational");
936 x2 = NumberMalloc(
"Norm-MakeRational");
937 if ( *arg2 == -SNUMBER ) {
938 if ( arg2[1] <= 1 ) {
939 NumberFree(x2,
"Norm-MakeRational");
940 NumberFree(x1,
"Norm-MakeRational");
943 else { x2[0] = arg2[1]; nx2 = 1; }
945 else if ( *arg2 > 0 ) {
947 nx2 = (ttstop[-1]-1)/2;
948 tc = arg2+ARGHEAD+1+nx2;
950 NumberFree(x2,
"Norm-MakeRational");
951 NumberFree(x1,
"Norm-MakeRational");
954 for ( i = 1; i < nx2; i++ )
if ( tc[i] != 0 ) {
955 NumberFree(x2,
"Norm-MakeRational");
956 NumberFree(x1,
"Norm-MakeRational");
960 for ( i = 0; i < nx2; i++ ) x2[i] = tc[i];
963 NumberFree(x2,
"Norm-MakeRational");
964 NumberFree(x1,
"Norm-MakeRational");
967 if ( BigLong(x1,ABS(nx1),x2,nx2) >= 0 ) {
968 UWORD *x3 = NumberMalloc(
"Norm-MakeRational");
969 UWORD *x4 = NumberMalloc(
"Norm-MakeRational");
971 DivLong(x1,nx1,x2,nx2,x3,&nx3,x4,&nx4);
972 for ( i = 0; i < ABS(nx4); i++ ) x1[i] = x4[i];
974 NumberFree(x4,
"Norm-MakeRational");
975 NumberFree(x3,
"Norm-MakeRational");
977 xx = (UWORD *)(TermMalloc(
"Norm-MakeRational"));
978 if ( MakeLongRational(BHEAD x1,nx1,x2,nx2,xx,&nxx) ) {
979 static int warnflag = 1;
981 MesPrint(
"%w Warning: fraction could not be reconstructed in MakeRational_");
984 ncoef = REDLENG(ncoef);
985 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,x1,nx1) )
989 ncoef = REDLENG(ncoef);
990 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,xx,nxx,
991 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
993 ncoef = INCLENG(ncoef);
994 TermFree(xx,
"Norm-MakeRational");
995 NumberFree(x2,
"Norm-MakeRational");
996 NumberFree(x1,
"Norm-MakeRational");
1000 if ( t[1] == FUNHEAD && AN.cTerm ) {
1001 ANsr = r; ANsm = m; ANsc = AN.cTerm;
1005 ncoef = REDLENG(ncoef);
1006 nnum = REDLENG(m[-1]);
1008 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,nnum,
1009 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
1010 ncoef = INCLENG(ncoef);
1015 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1016 if ( GetFirstBracket(termout,t[FUNHEAD+1]) < 0 )
goto FromNorm;
1017 if ( *termout == 0 )
goto NormZero;
1018 if ( *termout > 4 ) {
1020 while ( r < m ) *t++ = *r++;
1022 r2 = termout + *termout; r2 -= ABS(r2[-1]);
1023 while ( r < r1 ) *r2++ = *r++;
1025 while ( r3 < r2 ) *t++ = *r3++; *term = t - term;
1026 if ( AT.WorkPointer > term && AT.WorkPointer < t )
1034 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1037 POSITION oldondisk = AS.OldOnFile[t[FUNHEAD+1]];
1038 if ( e->replace == NEWLYDEFINEDEXPRESSION ) {
1039 AS.OldOnFile[t[FUNHEAD+1]] = e->onfile;
1041 if ( *t == FIRSTTERM ) {
1042 if ( GetFirstTerm(termout,t[FUNHEAD+1]) < 0 )
goto FromNorm;
1044 else if ( *t == CONTENTTERM ) {
1045 if ( GetContent(termout,t[FUNHEAD+1]) < 0 )
goto FromNorm;
1047 AS.OldOnFile[t[FUNHEAD+1]] = oldondisk;
1048 if ( *termout == 0 )
goto NormZero;
1052 WORD *r1, *r2, *r3, *r4, *r5, nr1, *rterm;
1053 r2 = termout + *termout; lnum = r2 - ABS(r2[-1]);
1054 nnum = REDLENG(r2[-1]);
1056 r1 = term + *term; r3 = r1 - ABS(r1[-1]);
1057 nr1 = REDLENG(r1[-1]);
1058 if ( Mully(BHEAD (UWORD *)lnum,&nnum,(UWORD *)r3,nr1) )
goto FromNorm;
1059 nnum = INCLENG(nnum); nr1 = ABS(nnum); lnum[nr1-1] = nnum;
1060 rterm = TermMalloc(
"FirstTerm/ContentTerm");
1061 r4 = rterm+1; r5 = term+1;
while ( r5 < t ) *r4++ = *r5++;
1062 r5 = termout+1;
while ( r5 < lnum ) *r4++ = *r5++;
1063 r5 = r;
while ( r5 < r3 ) *r4++ = *r5++;
1064 r5 = lnum; NCOPY(r4,r5,nr1);
1066 nr1 = *rterm; r1 = term; r2 = rterm; NCOPY(r1,r2,nr1);
1067 TermFree(rterm,
"FirstTerm/ContentTerm");
1068 if ( AT.WorkPointer > term && AT.WorkPointer < r1 )
1069 AT.WorkPointer = r1;
1073 else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
1074 DOLLARS d = Dollars + t[FUNHEAD+1], newd = 0;
1077 int nummodopt, dtype = -1;
1078 if ( AS.MultiThreaded && ( AC.mparallelflag == PARALLELFLAG ) ) {
1079 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
1080 if ( t[FUNHEAD+1] == ModOptdollars[nummodopt].number )
break;
1082 if ( nummodopt < NumModOptdollars ) {
1083 dtype = ModOptdollars[nummodopt].type;
1084 if ( dtype == MODLOCAL ) {
1085 d = ModOptdollars[nummodopt].dstruct+AT.identity;
1090 if ( d->where && ( d->type == DOLTERMS || d->type == DOLNUMBER ) ) {
1094 if ( ( newd = DolToTerms(BHEAD t[FUNHEAD+1]) ) == 0 )
1097 if ( newd->where[0] == 0 ) {
1098 M_free(newd,
"Copy of dollar variable");
1101 if ( *t == FIRSTTERM ) {
1102 idol = newd->where[0];
1103 for ( ido = 0; ido < idol; ido++ ) termout[ido] = newd->where[ido];
1105 else if ( *t == CONTENTTERM ) {
1107 tterm = newd->where;
1109 for ( ido = 0; ido < idol; ido++ ) termout[ido] = tterm[ido];
1112 if ( ContentMerge(BHEAD termout,tterm) < 0 )
goto FromNorm;
1117 if ( newd->factors ) M_free(newd->factors,
"Dollar factors");
1118 M_free(newd,
"Copy of dollar variable");
1126 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1127 x = TermsInExpression(t[FUNHEAD+1]);
1128 multermnum:
if ( x == 0 )
goto NormZero;
1131 if ( x > (LONG)WORDMASK ) { lnum[0] = x & WORDMASK;
1132 lnum[1] = x >> BITSINWORD; nnum = -2;
1134 else { lnum[0] = x; nnum = -1; }
1136 else if ( x > (LONG)WORDMASK ) {
1137 lnum[0] = x & WORDMASK;
1138 lnum[1] = x >> BITSINWORD;
1141 else { lnum[0] = x; nnum = 1; }
1142 ncoef = REDLENG(ncoef);
1143 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1145 ncoef = INCLENG(ncoef);
1147 else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
1148 x = TermsInDollar(t[FUNHEAD+1]);
1151 else { pcom[ncom++] = t; }
1155 case PATTERNFUNCTION:
1162 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1163 && t[FUNHEAD+1] >= 0 && t[FUNHEAD+2] == -SNUMBER
1164 && t[FUNHEAD+3] >= 0 && t[FUNHEAD+1] >= t[FUNHEAD+3] ) {
1165 if ( t[FUNHEAD+1] > t[FUNHEAD+3] ) {
1166 if ( GetBinom((UWORD *)lnum,&nnum,
1167 t[FUNHEAD+1],t[FUNHEAD+3]) )
goto FromNorm;
1168 if ( nnum == 0 )
goto NormZero;
1172 else pcom[ncom++] = t;
1178 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1179 if ( ( t[FUNHEAD+1] & 1 ) != 0 ) ncoef = -ncoef;
1181 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] )
1182 && ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) ) {
1183 UWORD *numer1,*denom1;
1184 WORD nsize = abs(t[t[1]-1]), nnsize, isize;
1185 nnsize = (nsize-1)/2;
1186 numer1 = (UWORD *)(t + FUNHEAD+ARGHEAD+1);
1187 denom1 = numer1 + nnsize;
1188 for ( isize = 1; isize < nnsize; isize++ ) {
1189 if ( denom1[isize] )
break;
1191 if ( ( denom1[0] != 1 ) || isize < nnsize ) {
1195 if ( ( numer1[0] & 1 ) != 0 ) ncoef = -ncoef;
1209 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1210 if ( t[FUNHEAD+1] < 0 ) ncoef = -ncoef;
1212 else if ( ( t[1] == FUNHEAD+2 ) && ( t[FUNHEAD] == -SYMBOL )
1213 && ( ( t[FUNHEAD+1] <= NumSymbols && t[FUNHEAD+1] > -MAXPOWER )
1214 && ( symbols[t[FUNHEAD+1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) ) {
1224 *m %= symbols[k].maxpower;
1225 if ( ( symbols[k].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
1226 if ( ( ( symbols[k].maxpower & 1 ) == 0 ) &&
1227 ( *m >= symbols[k].maxpower/2 ) ) {
1228 *m -= symbols[k].maxpower/2; ncoef = -ncoef;
1234 { *m = m[2]; m++; *m = m[2]; m++; i++; }
1240 }
while ( k < *m && --i > 0 );
1243 { m--; m[2] = *m; m--; m[2] = *m; i++; }
1251 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] ) ) {
1252 if ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) {
1253 if ( t[t[1]-1] < 0 ) ncoef = -ncoef;
1258 else if ( ( t[FUNHEAD+ARGHEAD]+FUNHEAD+ARGHEAD == t[1] )
1259 && ( t[FUNHEAD+ARGHEAD+1] == SYMBOL ) ) {
1260 WORD *ts = t + FUNHEAD+ARGHEAD+3;
1261 WORD its = ts[-1]-2;
1263 if ( ( *ts != 0 ) && ( ( *ts > NumSymbols || *ts <= -MAXPOWER )
1264 || ( symbols[*ts].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY ) ) {
1273 if ( t[t[1]-1] < 0 ) ncoef = -ncoef;
1274 ts = t + FUNHEAD+ARGHEAD+3;
1285 if ( ( ts[-1] <= NumSymbols && ts[-1] > -MAXPOWER ) &&
1286 ( symbols[ts[-1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
1287 *m %= symbols[ts[-1]].maxpower;
1288 if ( *m < 0 ) *m += symbols[ts[-1]].maxpower;
1289 if ( ( symbols[ts[-1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
1290 if ( ( ( symbols[ts[-1]].maxpower & 1 ) == 0 ) &&
1291 ( *m >= symbols[ts[-1]].maxpower/2 ) ) {
1292 *m -= symbols[ts[-1]].maxpower/2; ncoef = -ncoef;
1299 { *m = m[2]; m++; *m = m[2]; m++; i++; }
1306 }
while ( *ts < *m && --i > 0 );
1309 { m--; m[2] = *m; m--; m[2] = *m; i++; }
1319 signogood: pcom[ncom++] = t;
1322 else pcom[ncom++] = t;
1329 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1331 if ( k < 0 ) k = -k;
1332 if ( k == 0 )
goto NormZero;
1333 *((UWORD *)lnum) = k; nnum = 1;
1337 else if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SYMBOL ) {
1339 if ( ( k > NumSymbols || k <= -MAXPOWER )
1340 || ( symbols[k].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
1343 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] )
1344 && ( t[1] == FUNHEAD+ARGHEAD+t[FUNHEAD+ARGHEAD] ) ) {
1345 if ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) {
1347 absnosymbols: ts = t + t[1] -1;
1348 ncoef = REDLENG(ncoef);
1349 nnum = REDLENG(*ts);
1350 if ( nnum < 0 ) nnum = -nnum;
1351 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,
1352 (UWORD *)(ts-ABS(*ts)+1),nnum,
1353 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
1354 ncoef = INCLENG(ncoef);
1359 else if ( t[FUNHEAD+ARGHEAD+1] == SYMBOL ) {
1360 WORD *ts = t+FUNHEAD+ARGHEAD+1;
1361 WORD its = ts[1] - 2;
1365 else if ( ( *ts > NumSymbols || *ts <= -MAXPOWER )
1366 || ( symbols[*ts].complex & VARTYPEROOTOFUNITY )
1367 != VARTYPEROOTOFUNITY )
goto absnogood;
1373 absnogood: pcom[ncom++] = t;
1376 else pcom[ncom++] = t;
1384 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1385 && t[FUNHEAD+2] == -SNUMBER && t[FUNHEAD+3] > 1 ) {
1387 tmod = (t[FUNHEAD+1]%t[FUNHEAD+3]);
1388 if ( tmod < 0 ) tmod += t[FUNHEAD+3];
1389 if ( *t == MOD2FUNCTION && tmod > t[FUNHEAD+3]/2 )
1390 tmod -= t[FUNHEAD+3];
1392 *((UWORD *)lnum) = -tmod;
1395 else if ( tmod > 0 ) {
1396 *((UWORD *)lnum) = tmod;
1402 else if ( t[1] > t[FUNHEAD+2] && t[FUNHEAD] > 0
1403 && t[FUNHEAD+t[FUNHEAD]] == -SNUMBER
1404 && t[FUNHEAD+t[FUNHEAD]+1] > 1
1405 && t[1] == FUNHEAD+2+t[FUNHEAD] ) {
1406 WORD *ttt = t+FUNHEAD, iii;
1408 if ( *ttt == ttt[ARGHEAD]+ARGHEAD &&
1409 ttt[ARGHEAD] == ABS(iii)+1 ) {
1411 WORD cmod = ttt[*ttt+1];
1413 if ( *t == MODFUNCTION ) {
1414 if ( TakeModulus((UWORD *)(ttt+ARGHEAD+1)
1415 ,&iii,(UWORD *)(&cmod),ncmod,UNPACK|NOINVERSES) )
1419 if ( TakeModulus((UWORD *)(ttt+ARGHEAD+1)
1420 ,&iii,(UWORD *)(&cmod),ncmod,UNPACK|POSNEG|NOINVERSES) )
1423 if ( *t == MOD2FUNCTION && ttt[ARGHEAD+1] > cmod/2 && iii > 0 ) {
1424 ttt[ARGHEAD+1] -= cmod;
1426 if ( ttt[ARGHEAD+1] < 0 ) {
1427 *((UWORD *)lnum) = -ttt[ARGHEAD+1];
1430 else if ( ttt[ARGHEAD+1] > 0 ) {
1431 *((UWORD *)lnum) = ttt[ARGHEAD+1];
1438 else if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1439 *((UWORD *)lnum) = t[FUNHEAD+1];
1440 if ( *lnum == 0 )
goto NormZero;
1444 else if ( ( ( t[FUNHEAD] < 0 ) && ( t[FUNHEAD] == -SNUMBER )
1445 && ( t[1] >= ( FUNHEAD+6+ARGHEAD ) )
1446 && ( t[FUNHEAD+2] >= 4+ARGHEAD )
1447 && ( t[t[1]-1] == t[FUNHEAD+2+ARGHEAD]-1 ) ) ||
1448 ( ( t[FUNHEAD] > 0 )
1449 && ( t[FUNHEAD]-ARGHEAD-1 == ABS(t[FUNHEAD+t[FUNHEAD]-1]) )
1450 && ( t[FUNHEAD+t[FUNHEAD]]-ARGHEAD-1 == t[t[1]-1] ) ) ) {
1454 WORD *ttt = t + t[1], iii, iii1;
1455 UWORD coefbuf[2], *coef2, ncoef2;
1456 iii = (ttt[-1]-1)/2;
1458 if ( ttt[-1] != 1 ) {
1464 for ( iii1 = 0; iii1 < iii; iii1++ ) {
1465 if ( ttt[iii1] != 0 )
goto exitfromhere;
1476 nnum = -1; lnum[0] = -ttt[1]; lnum[1] = 1;
1479 nnum = 1; lnum[0] = ttt[1]; lnum[1] = 1;
1483 nnum = ABS(ttt[ttt[0]-1] - 1);
1484 for ( iii = 0; iii < nnum; iii++ ) {
1485 lnum[iii] = ttt[ARGHEAD+1+iii];
1488 if ( ttt[ttt[0]-1] < 0 ) nnum = -nnum;
1493 ncoef2 = 3; *coef2 = (UWORD)(ttt[1]);
1497 coef2 = (UWORD *)(ttt+ARGHEAD+1);
1498 ncoef2 = (ttt[ttt[0]-1]-1)/2;
1500 if ( TakeModulus((UWORD *)lnum,&nnum,(UWORD *)coef2,ncoef2,
1501 UNPACK|NOINVERSES|FROMFUNCTION) ) {
1504 if ( *t == MOD2FUNCTION && nnum > 0 ) {
1505 UWORD *coef3 = NumberMalloc(
"Mod2Function"), two = 2;
1507 if ( MulLong((UWORD *)lnum,nnum,&two,1,coef3,&ncoef3) )
1509 if ( BigLong(coef3,ncoef3,(UWORD *)coef2,ncoef2) > 0 ) {
1511 AddLong((UWORD *)lnum,nnum,(UWORD *)coef2,ncoef2
1512 ,(UWORD *)lnum,&nnum);
1515 NumberFree(coef3,
"Mod2Function");
1520 ncoef = REDLENG(ncoef);
1521 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1523 ncoef = INCLENG(ncoef);
1525 else pcom[ncom++] = t;
1529 WORD argcount = 0, *tc, *ts, xc, xs, *tcc;
1530 UWORD *Num1, *Num2, *Num3, *Num4;
1531 WORD size1, size2, size3, size4, space;
1532 tc = t+FUNHEAD; ts = t + t[1];
1533 while ( argcount < 3 && tc < ts ) { NEXTARG(tc); argcount++; }
1534 if ( argcount != 2 )
goto defaultcase;
1535 if ( t[FUNHEAD] == -SNUMBER ) {
1536 if ( t[FUNHEAD+1] <= 1 )
goto defaultcase;
1537 if ( t[FUNHEAD+2] == -SNUMBER ) {
1538 if ( t[FUNHEAD+3] <= 1 )
goto defaultcase;
1539 Num2 = NumberMalloc(
"modinverses");
1540 *Num2 = t[FUNHEAD+3]; size2 = 1;
1543 if ( ts[-1] < 0 )
goto defaultcase;
1544 if ( ts[-1] != t[FUNHEAD+2]-ARGHEAD-1 )
goto defaultcase;
1547 if ( *tcc != 1 )
goto defaultcase;
1548 for ( i = 1; i < xs; i++ ) {
1549 if ( tcc[i] != 0 )
goto defaultcase;
1551 Num2 = NumberMalloc(
"modinverses");
1553 for ( i = 0; i < xs; i++ ) Num2[i] = t[FUNHEAD+ARGHEAD+3+i];
1555 Num1 = NumberMalloc(
"modinverses");
1556 *Num1 = t[FUNHEAD+1]; size1 = 1;
1559 tc = t + FUNHEAD + t[FUNHEAD];
1560 if ( tc[-1] < 0 )
goto defaultcase;
1561 if ( tc[-1] != t[FUNHEAD]-ARGHEAD-1 )
goto defaultcase;
1564 if ( *tcc != 1 )
goto defaultcase;
1565 for ( i = 1; i < xc; i++ ) {
1566 if ( tcc[i] != 0 )
goto defaultcase;
1568 if ( *tc == -SNUMBER ) {
1569 if ( tc[1] <= 1 )
goto defaultcase;
1570 Num2 = NumberMalloc(
"modinverses");
1571 *Num2 = tc[1]; size2 = 1;
1574 if ( ts[-1] < 0 )
goto defaultcase;
1575 if ( ts[-1] != t[FUNHEAD+2]-ARGHEAD-1 )
goto defaultcase;
1578 if ( *tcc != 1 )
goto defaultcase;
1579 for ( i = 1; i < xs; i++ ) {
1580 if ( tcc[i] != 0 )
goto defaultcase;
1582 Num2 = NumberMalloc(
"modinverses");
1584 for ( i = 0; i < xs; i++ ) Num2[i] = tc[ARGHEAD+1+i];
1586 Num1 = NumberMalloc(
"modinverses");
1588 for ( i = 0; i < xc; i++ ) Num1[i] = t[FUNHEAD+ARGHEAD+1+i];
1590 Num3 = NumberMalloc(
"modinverses");
1591 Num4 = NumberMalloc(
"modinverses");
1592 GetLongModInverses(BHEAD Num1,size1,Num2,size2
1593 ,Num3,&size3,Num4,&size4);
1602 if ( ( size3 == 1 || size3 == -1 ) && (*Num3&TOPBITONLY) == 0 ) space += 2;
1603 else space += ARGHEAD + 2*ABS(size3) + 2;
1604 if ( ( size4 == 1 || size4 == -1 ) && (*Num4&TOPBITONLY) == 0 ) space += 2;
1605 else space += ARGHEAD + 2*ABS(size4) + 2;
1606 tt = term + *term; u = tt + space;
1607 while ( tt >= ts ) *--u = *--tt;
1608 m += space; r += space;
1611 if ( ( size3 == 1 || size3 == -1 ) && (*Num3&TOPBITONLY) == 0 ) {
1612 *ts++ = -SNUMBER; *ts = (WORD)(*Num3);
1613 if ( size3 < 0 ) *ts = -*ts;
1617 *ts++ = 2*ABS(size3)+ARGHEAD+2;
1618 *ts++ = 0; FILLARG(ts)
1619 *ts++ = 2*ABS(size3)+1;
1620 for ( i = 0; i < ABS(size3); i++ ) *ts++ = Num3[i];
1622 for ( i = 1; i < ABS(size3); i++ ) *ts++ = 0;
1623 if ( size3 < 0 ) *ts++ = 2*size3-1;
1624 else *ts++ = 2*size3+1;
1626 if ( ( size4 == 1 || size4 == -1 ) && (*Num4&TOPBITONLY) == 0 ) {
1627 *ts++ = -SNUMBER; *ts = *Num4;
1628 if ( size4 < 0 ) *ts = -*ts;
1632 *ts++ = 2*ABS(size4)+ARGHEAD+2;
1633 *ts++ = 0; FILLARG(ts)
1634 *ts++ = 2*ABS(size4)+2;
1635 for ( i = 0; i < ABS(size4); i++ ) *ts++ = Num4[i];
1637 for ( i = 1; i < ABS(size4); i++ ) *ts++ = 0;
1638 if ( size4 < 0 ) *ts++ = 2*size4-1;
1639 else *ts++ = 2*size4+1;
1641 NumberFree(Num4,
"modinverses");
1642 NumberFree(Num3,
"modinverses");
1643 NumberFree(Num1,
"modinverses");
1644 NumberFree(Num2,
"modinverses");
1651 #ifdef NEWGCDFUNCTION 1657 WORD *num1, *num2, size1, size2, stor1, stor2, *ttt, ti;
1658 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1659 && t[FUNHEAD+2] == -SNUMBER && t[FUNHEAD+1] != 0
1660 && t[FUNHEAD+3] != 0 ) {
1661 stor1 = t[FUNHEAD+1];
1662 stor2 = t[FUNHEAD+3];
1663 if ( stor1 < 0 ) stor1 = -stor1;
1664 if ( stor2 < 0 ) stor2 = -stor2;
1665 num1 = &stor1; num2 = &stor2;
1669 else if ( t[1] > FUNHEAD+4 ) {
1670 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] != 0
1671 && t[FUNHEAD+2] == t[1]-FUNHEAD-2 &&
1672 ABS(t[t[1]-1]) == t[FUNHEAD+2]-1-ARGHEAD ) {
1674 size2 = ABS(num2[-1]);
1677 size2 = (size2-1)/2;
1679 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1680 if ( ti == 1 && ttt[-1] == 1 ) {
1681 stor1 = t[FUNHEAD+1];
1682 if ( stor1 < 0 ) stor1 = -stor1;
1687 else pcom[ncom++] = t;
1689 else if ( t[FUNHEAD] > 0 &&
1690 t[FUNHEAD]-1-ARGHEAD == ABS(t[t[FUNHEAD]+FUNHEAD-1]) ) {
1691 num1 = t + FUNHEAD + t[FUNHEAD];
1692 size1 = ABS(num1[-1]);
1695 size1 = (size1-1)/2;
1697 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1698 if ( ti == 1 && ttt[-1] == 1 ) {
1699 if ( t[1]-FUNHEAD == t[FUNHEAD]+2 && t[t[1]-2] == -SNUMBER
1700 && t[t[1]-1] != 0 ) {
1702 if ( stor2 < 0 ) stor2 = -stor2;
1707 else if ( t[1]-FUNHEAD == t[FUNHEAD]+t[FUNHEAD+t[FUNHEAD]]
1708 && ABS(t[t[1]-1]) == t[FUNHEAD+t[FUNHEAD]] - ARGHEAD-1 ) {
1710 size2 = ABS(num2[-1]);
1713 size2 = (size2-1)/2;
1715 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1716 if ( ti == 1 && ttt[-1] == 1 ) {
1717 gcdcalc:
if ( GcdLong(BHEAD (UWORD *)num1,size1,(UWORD *)num2,size2
1718 ,(UWORD *)lnum,&nnum) )
goto FromNorm;
1721 else pcom[ncom++] = t;
1723 else pcom[ncom++] = t;
1725 else pcom[ncom++] = t;
1727 else pcom[ncom++] = t;
1729 else pcom[ncom++] = t;
1733 WORD *gcd = AT.WorkPointer;
1734 if ( ( gcd = EvaluateGcd(BHEAD t) ) == 0 )
goto FromNorm;
1735 if ( *gcd == 4 && gcd[1] == 1 && gcd[2] == 1 && gcd[4] == 0 ) {
1736 AT.WorkPointer = gcd;
1738 else if ( gcd[*gcd] == 0 ) {
1739 WORD *t1, iii, change, *num, *den, numsize, densize;
1740 if ( gcd[*gcd-1] < *gcd-1 ) {
1742 for ( iii = 2; iii < t1[1]; iii += 2 ) {
1743 change = ExtraSymbol(t1[iii],t1[iii+1],nsym,ppsym,&ncoef);
1745 ppsym += change << 1;
1749 iii = t1[-1]; num = t1-iii; numsize = (iii-1)/2;
1750 den = num + numsize; densize = numsize;
1751 while ( numsize > 1 && num[numsize-1] == 0 ) numsize--;
1752 while ( densize > 1 && den[densize-1] == 0 ) densize--;
1753 if ( numsize > 1 || num[0] != 1 ) {
1754 ncoef = REDLENG(ncoef);
1755 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)num,numsize) )
goto FromNorm;
1756 ncoef = INCLENG(ncoef);
1758 if ( densize > 1 || den[0] != 1 ) {
1759 ncoef = REDLENG(ncoef);
1760 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)den,densize) )
goto FromNorm;
1761 ncoef = INCLENG(ncoef);
1763 AT.WorkPointer = gcd;
1777 LONG size = AT.WorkPointer - gcd;
1779 ss =
AddRHS(AT.ebufnum,1);
1783 C->
rhs[C->numrhs+1] = ss;
1786 t[0] = SUBEXPRESSION;
1793 while ( r < tt ) *t++ = *r++;
1802 MesPrint(
" Unexpected call to EvaluateGCD");
1808 if ( t[1] == FUNHEAD )
break;
1810 WORD *ttt = t + FUNHEAD;
1811 WORD *tttstop = t + t[1];
1813 while ( ttt < tttstop ) {
1815 if ( ttt[ARGHEAD]-1 > ABS(ttt[*ttt-1]) )
goto nospec;
1819 if ( *ttt != -SNUMBER )
goto nospec;
1830 for ( iii = 0; iii < ttt[ARGHEAD]; iii++ )
1831 n_llnum[iii] = ttt[ARGHEAD+iii];
1836 if ( ttt[1] == 0 ) {
1837 n_llnum[0] = n_llnum[1] = n_llnum[2] = n_llnum[3] = 0;
1841 if ( ttt[1] > 0 ) { n_llnum[1] = ttt[1]; n_llnum[3] = 3; }
1842 else { n_llnum[1] = -ttt[1]; n_llnum[3] = -3; }
1850 while ( ttt < tttstop ) {
1852 if ( n_llnum[0] == 0 ) {
1853 if ( ( *t == MINFUNCTION && ttt[*ttt-1] < 0 )
1854 || ( *t == MAXFUNCTION && ttt[*ttt-1] > 0 ) )
1860 if ( ( iii > 0 && *t == MINFUNCTION )
1861 || ( iii < 0 && *t == MAXFUNCTION ) ) {
1862 for ( iii = 0; iii < ttt[0]; iii++ )
1863 n_llnum[iii] = ttt[iii];
1869 if ( n_llnum[0] == 0 ) {
1870 if ( ( *t == MINFUNCTION && ttt[1] < 0 )
1871 || ( *t == MAXFUNCTION && ttt[1] > 0 ) )
1874 else if ( ttt[1] == 0 ) {
1875 if ( ( *t == MINFUNCTION && n_llnum[*n_llnum-1] > 0 )
1876 || ( *t == MAXFUNCTION && n_llnum[*n_llnum-1] < 0 ) ) {
1881 tterm[0] = 4; tterm[2] = 1;
1882 if ( ttt[1] < 0 ) { tterm[1] = -ttt[1]; tterm[3] = -3; }
1883 else { tterm[1] = ttt[1]; tterm[3] = 3; }
1885 if ( ( iii > 0 && *t == MINFUNCTION )
1886 || ( iii < 0 && *t == MAXFUNCTION ) ) {
1887 for ( iii = 0; iii < 4; iii++ )
1888 n_llnum[iii] = tterm[iii];
1894 if ( n_llnum[0] == 0 )
goto NormZero;
1895 ncoef = REDLENG(ncoef);
1896 nnum = REDLENG(n_llnum[*n_llnum-1]);
1897 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)lnum,nnum,
1898 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
1899 ncoef = INCLENG(ncoef);
1902 case INVERSEFACTORIAL:
1903 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] >= 0 ) {
1904 if ( Factorial(BHEAD t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
1906 ncoef = REDLENG(ncoef);
1907 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
goto FromNorm;
1908 ncoef = INCLENG(ncoef);
1911 nospec: pcom[ncom++] = t;
1915 if ( ( t[FUNHEAD] == -SYMBOL )
1916 && ( t[FUNHEAD+1] > 0 ) && ( t[1] == FUNHEAD+2 ) ) {
1917 *((UWORD *)lnum) = symbols[t[FUNHEAD+1]].maxpower;
1921 else { pcom[ncom++] = t; }
1924 if ( ( t[FUNHEAD] == -SYMBOL )
1925 && ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+2 ) ) {
1926 *((UWORD *)lnum) = symbols[t[FUNHEAD+1]].minpower;
1930 else { pcom[ncom++] = t; }
1933 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER
1934 && t[FUNHEAD+1] > 0 ) {
1935 UWORD xp = (UWORD)(
NextPrime(BHEAD t[FUNHEAD+1]));
1936 ncoef = REDLENG(ncoef);
1937 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,&xp,1) )
goto FromNorm;
1938 ncoef = INCLENG(ncoef);
1940 else goto defaultcase;
1943 ncoef = REDLENG(ncoef);
1944 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(t+3),t[2]) )
goto FromNorm;
1945 ncoef = INCLENG(ncoef);
1950 if ( t[3] & 1 ) ncoef = -ncoef;
1952 else if ( t[2] == 0 ) {
1953 if ( t[3] < 0 )
goto NormInf;
1958 if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) )
goto FromNorm;
1959 ncoef = REDLENG(ncoef);
1961 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1964 else if ( t[3] > 0 ) {
1965 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1968 ncoef = INCLENG(ncoef);
1975 if ( t[1] == FUNHEAD ) {
1976 MLOCK(ErrorMessageLock);
1977 MesPrint(
"Gamma matrix without spin line encountered.");
1978 MUNLOCK(ErrorMessageLock);
1986 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG ) {
1988 t[2] |= DIRTYSYMFLAG;
1991 ScanCont:
while ( t < r ) {
1992 if ( *t >= AM.OffsetIndex &&
1993 ( *t >= AM.DumInd || ( *t < AM.WilInd &&
1994 indices[*t-AM.OffsetIndex].dimension ) ) )
2004 if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) )
2006 if ( *rr == -SNUMBER && rr[1] == 1 )
break;
2007 if ( *rr <= -FUNCTION ) k = *rr;
2009 if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) ) {
2010 if ( k == 0 )
goto NormZZ;
2013 if ( *rr == -SNUMBER && rr[1] > 0 && rr[1] < MAXPOWER && k < 0 ) {
2015 if ( functions[k-FUNCTION].commute ) {
2016 for ( i = 0; i < rr[1]; i++ ) pnco[nnco++] = rr-1;
2019 for ( i = 0; i < rr[1]; i++ ) pcom[ncom++] = rr-1;
2023 if ( k == 0 )
goto NormZero;
2024 if ( t[FUNHEAD] == -SYMBOL && *rr == -SNUMBER && t[1] == FUNHEAD+4 ) {
2025 if ( rr[1] < MAXPOWER ) {
2026 t[FUNHEAD+2] = t[FUNHEAD+1]; t += FUNHEAD+2;
2036 t[2] &= ~DIRTYSYMFLAG;
2042 t[2] &= ~DIRTYSYMFLAG;
2049 if ( *t == 0 )
goto NormZero;
2050 if ( *t > 0 && *t < AM.OffsetIndex ) {
2053 ncoef = REDLENG(ncoef);
2054 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2056 ncoef = INCLENG(ncoef);
2058 else if ( *t == NOINDEX ) t++;
2059 else pind[nind++] = *t++;
2062 case SUBEXPRESSION :
2063 if ( t[3] == 0 )
break;
2075 if ( t[2] == 0 )
goto defaultcase;
2076 if ( t[FUNHEAD] != -SNUMBER || t[FUNHEAD+1] < 0 )
goto defaultcase;
2077 if ( t[FUNHEAD+2] == -SNUMBER ) {
2078 if ( t[FUNHEAD+1] == 0 && t[FUNHEAD+3] == 0 )
goto NormZZ;
2079 if ( t[FUNHEAD+1] == 0 )
break;
2080 if ( t[FUNHEAD+3] < 0 ) {
2081 AT.WorkPointer[0] = -t[FUNHEAD+3];
2085 AT.WorkPointer[0] = t[FUNHEAD+3];
2088 AT.WorkPointer[1] = 1;
2090 else if ( t[FUNHEAD+2] == t[1]-FUNHEAD-2
2091 && t[FUNHEAD+2] == t[FUNHEAD+2+ARGHEAD]+ARGHEAD
2092 && ABS(t[t[1]-1]) == t[FUNHEAD+2+ARGHEAD] - 1 ) {
2094 if ( t[FUNHEAD+1] == 0 )
break;
2095 i = t[t[1]-1]; r1 = t + FUNHEAD+ARGHEAD+3;
2098 r2 = AT.WorkPointer;
2099 while ( --i >= 0 ) *r2++ = *r1++;
2101 else goto defaultcase;
2102 if ( TakeRatRoot((UWORD *)AT.WorkPointer,&nc,t[FUNHEAD+1]) ) {
2106 ncoef = REDLENG(ncoef);
2107 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,nc) )
2109 if ( nc < 0 ) nc = -nc;
2110 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(AT.WorkPointer+nc),nc) )
2112 ncoef = INCLENG(ncoef);
2115 case RANDOMFUNCTION :
2117 WORD nnc, nc, nca, nr;
2128 if ( t[1] == FUNHEAD )
goto defaultcase;
2129 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER &&
2130 t[FUNHEAD+1] > 0 ) {
2131 if ( t[FUNHEAD+1] == 1 )
break;
2133 ((UWORD *)AT.WorkPointer)[0] = wranf(BHEAD0);
2134 ((UWORD *)AT.WorkPointer)[1] = wranf(BHEAD0);
2136 if ( ((UWORD *)AT.WorkPointer)[1] == 0 ) {
2138 if ( ((UWORD *)AT.WorkPointer)[0] == 0 ) {
2142 xx = (UWORD)(t[FUNHEAD+1]);
2144 DivLong((UWORD *)AT.WorkPointer,nr
2146 ,((UWORD *)AT.WorkPointer)+4,&nnc
2147 ,((UWORD *)AT.WorkPointer)+2,&nc);
2148 ((UWORD *)AT.WorkPointer)[4] = 0;
2149 ((UWORD *)AT.WorkPointer)[5] = 0;
2150 ((UWORD *)AT.WorkPointer)[6] = 1;
2151 DivLong((UWORD *)AT.WorkPointer+4,3
2153 ,((UWORD *)AT.WorkPointer)+9,&nnc
2154 ,((UWORD *)AT.WorkPointer)+7,&nca);
2155 AddLong((UWORD *)AT.WorkPointer+4,3
2156 ,((UWORD *)AT.WorkPointer)+7,-nca
2157 ,((UWORD *)AT.WorkPointer)+9,&nnc);
2158 if ( BigLong((UWORD *)AT.WorkPointer,nr
2159 ,((UWORD *)AT.WorkPointer)+9,nnc) >= 0 )
goto redoshort;
2163 AT.WorkPointer[2] = (WORD)xx;
2166 ncoef = REDLENG(ncoef);
2167 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+2,nc) )
2169 ncoef = INCLENG(ncoef);
2171 else if ( t[FUNHEAD] > 0 && t[1] == t[FUNHEAD]+FUNHEAD
2172 && ABS(t[t[1]-1]) == t[FUNHEAD]-1-ARGHEAD && t[t[1]-1] > 0 ) {
2173 WORD nna, nnb, nni, nnb2, nnb2a;
2178 nnt = (UWORD *)(t+t[1]-1-nnb);
2179 if ( *nnt != 1 )
goto defaultcase;
2180 for ( nni = 1; nni < nnb; nni++ ) {
2181 if ( nnt[nni] != 0 )
goto defaultcase;
2183 nnt = (UWORD *)(t + FUNHEAD + ARGHEAD + 1);
2185 for ( nni = 0; nni < nnb2; nni++ ) {
2186 ((UWORD *)AT.WorkPointer)[nni] = wranf(BHEAD0);
2189 while ( nnb2a > 0 && ((UWORD *)AT.WorkPointer)[nnb2a-1] == 0 ) nnb2a--;
2191 DivLong((UWORD *)AT.WorkPointer,nnb2a
2193 ,((UWORD *)AT.WorkPointer)+2*nnb2,&nnc
2194 ,((UWORD *)AT.WorkPointer)+nnb2,&nc);
2195 for ( nni = 0; nni < nnb2; nni++ ) {
2196 ((UWORD *)AT.WorkPointer)[nni+2*nnb2] = 0;
2198 ((UWORD *)AT.WorkPointer)[3*nnb2] = 1;
2199 DivLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2201 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc
2202 ,((UWORD *)AT.WorkPointer)+3*nnb2+1,&nca);
2203 AddLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2204 ,((UWORD *)AT.WorkPointer)+3*nnb2+1,-nca
2205 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc);
2206 if ( BigLong((UWORD *)AT.WorkPointer,nnb2a
2207 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,nnc) >= 0 )
goto redoshort;
2211 for ( nni = 0; nni < nnb; nni++ ) {
2212 ((UWORD *)AT.WorkPointer)[nnb2+nni] = nnt[nni];
2216 ncoef = REDLENG(ncoef);
2217 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+nnb2,nc) )
2219 ncoef = INCLENG(ncoef);
2221 else goto defaultcase;
2225 if ( *t == RANPERM && t[1] > FUNHEAD && t[FUNHEAD] <= -FUNCTION ) {
2227 WORD *mm, *ww, *ow = AT.WorkPointer;
2228 WORD *Array, *targ, *argstop, narg = 0, itot;
2232 while ( targ < argstop ) {
2233 narg++; NEXTARG(targ);
2235 WantAddPointers(narg);
2236 pwork = AT.pWorkSpace + AT.pWorkPointer;
2237 targ = t+FUNHEAD+1; narg = 0;
2238 while ( targ < argstop ) {
2239 pwork[narg++] = targ;
2246 ow = AT.WorkPointer;
2247 Array = AT.WorkPointer;
2248 AT.WorkPointer += narg;
2249 for ( i = 0; i < narg; i++ ) Array[i] = i;
2250 for ( i = 2; i <= narg; i++ ) {
2251 itot = (WORD)(iranf(BHEAD i));
2252 for ( j = 0; j < itot; j++ ) CYCLE1(WORD,Array,i)
2254 mm = AT.WorkPointer;
2255 *mm++ = -t[FUNHEAD];
2257 for ( ie = 2; ie < FUNHEAD; ie++ ) *mm++ = t[ie];
2258 for ( i = 0; i < narg; i++ ) {
2259 ww = pwork[Array[i]];
2262 mm = AT.WorkPointer; t++; ww = t;
2263 i = mm[1]; NCOPY(ww,mm,i)
2264 AT.WorkPointer = ow;
2270 if ( ( t[2] & DIRTYFLAG ) != 0 && t[FUNHEAD] <= -FUNCTION
2271 && t[FUNHEAD+1] == -SNUMBER && t[FUNHEAD+2] > 0 ) {
2272 WORD *rr = t+t[1], *mm = t+FUNHEAD+3, *tt, *tt1, *tt2, num = 0;
2276 while ( mm < rr ) { num++; NEXTARG(mm); }
2277 if ( num < t[FUNHEAD+2] ) { pnco[nnco++] = t;
break; }
2278 *t = -t[FUNHEAD]; mm = t+FUNHEAD+3;
2280 while ( --i > 0 ) { NEXTARG(mm); }
2281 tt = TermMalloc(
"Select_");
2284 for ( i = 0; i < *mm; i++ ) *tt1++ = mm[i];
2286 else if ( *mm <= -FUNCTION ) { *tt1++ = *mm; }
2287 else { *tt1++ = mm[0]; *tt1++ = mm[1]; }
2289 while ( tt2 < mm ) *tt1++ = *tt2++;
2290 i = tt1-tt; tt1 = tt; tt2 = t+FUNHEAD;
2292 TermFree(tt,
"Select_");
2294 while ( mm < rr ) *tt2++ = *mm++;
2297 while ( mm < rr ) *tt2++ = *mm++;
2301 else pnco[nnco++] = t;
2312 if ( t[1] <= FUNHEAD )
break;
2317 if ( *to == ARGHEAD )
goto NormZero;
2321 if ( to[ARGHEAD] != j+1 )
goto NoInteg;
2322 if ( rr >= r ) k = -1;
2323 else if ( *rr == ARGHEAD ) { k = 0; rr += ARGHEAD; }
2324 else if ( *rr == -SNUMBER ) { k = rr[1]; rr += 2; }
2326 if ( rr != r )
goto NoInteg;
2327 if ( k > 1 || k < -1 )
goto NoInteg;
2330 i = ( i < 0 ) ? -j: j;
2331 UnPack((UWORD *)to,i,&den,&num);
2339 if ( AN.NoScrat2 == 0 ) {
2340 AN.NoScrat2 = (UWORD *)Malloc1((AM.MaxTal+2)*
sizeof(UWORD),
"Normalize");
2342 if ( DivLong((UWORD *)to,num,(UWORD *)(to+j),den
2343 ,(UWORD *)AT.WorkPointer,&num,AN.NoScrat2,&den) )
goto FromNorm;
2344 if ( k < 0 && den < 0 ) {
2347 if ( AddLong((UWORD *)AT.WorkPointer,num
2348 ,AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2351 else if ( k > 0 && den > 0 ) {
2354 if ( AddLong((UWORD *)AT.WorkPointer,num,
2355 AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2360 else if ( *to == -SNUMBER ) {
2361 if ( to[1] < 0 ) { *AT.WorkPointer = -to[1]; num = -1; }
2362 else if ( to[1] == 0 )
goto NormZero;
2363 else { *AT.WorkPointer = to[1]; num = 1; }
2366 if ( num == 0 )
goto NormZero;
2367 ncoef = REDLENG(ncoef);
2368 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,num) )
2370 ncoef = INCLENG(ncoef);
2379 if ( *t < FUNCTION ) {
2380 MLOCK(ErrorMessageLock);
2381 MesPrint(
"Illegal code in Norm");
2385 AO.OutFill = AO.OutputLine = OutBuf;
2390 while ( --i >= 0 ) {
2391 TalToLine((UWORD)(*t++));
2392 TokenToLine((UBYTE *)
" ");
2398 MUNLOCK(ErrorMessageLock);
2401 if ( *t == REPLACEMENT ) {
2402 if ( AR.Eside != LHSIDE ) ReplaceVeto--;
2410 if ( *t == DUMMYFUN || *t == DUMMYTEN ) {}
2412 if ( *t < (FUNCTION + WILDOFFSET) ) {
2413 if ( ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2414 || ( functions[*t-FUNCTION].minnumargs > 0 ) )
2415 && ( ( t[2] & DIRTYFLAG ) != 0 ) ) {
2419 WORD *ta = t + FUNHEAD, *tb = t + t[1];
2421 while ( ta < tb ) { numarg++; NEXTARG(ta) }
2422 if ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2423 && ( numarg >= functions[*t-FUNCTION].maxnumargs ) )
2425 if ( ( functions[*t-FUNCTION].minnumargs > 0 )
2426 && ( numarg < functions[*t-FUNCTION].minnumargs ) )
2430 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION].tabl == 0 ) ) {
2432 t[2] |= DIRTYSYMFLAG;
2434 if ( functions[*t-FUNCTION].commute ) { pnco[nnco++] = t; }
2435 else { pcom[ncom++] = t; }
2438 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION-WILDOFFSET].tabl == 0 ) ) {
2440 t[2] |= DIRTYSYMFLAG;
2442 if ( functions[*t-FUNCTION-WILDOFFSET].commute ) {
2445 else { pcom[ncom++] = t; }
2451 if ( ( *t < (FUNCTION + WILDOFFSET)
2452 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) || (
2453 *t >= (FUNCTION + WILDOFFSET)
2454 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION ) ) {
2455 if ( *t >= GAMMA && *t <= GAMMASEVEN ) t++;
2458 if ( *t >= AM.OffsetIndex && ( *t >= AM.DumInd
2459 || ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2462 else if ( *t == FUNNYWILD ) { t++; }
2477 else if ( *t <= -FUNCTION ) t++;
2478 else if ( *t == -INDEX ) {
2479 if ( t[1] >= AM.OffsetIndex &&
2480 ( t[1] >= AM.DumInd || ( t[1] < AM.WilInd
2481 && indices[t[1]-AM.OffsetIndex].dimension ) ) )
2485 else if ( *t == -SYMBOL ) {
2486 if ( t[1] >= MAXPOWER && t[1] < 2*MAXPOWER ) {
2490 else if ( t[1] < -MAXPOWER && t[1] > -2*MAXPOWER ) {
2506 r = t = ANsr; m = ANsm;
2507 ANsc = ANsm = ANsr = 0;
2520 for ( k = 0, i = 0; i < nden; i++ ) {
2522 if ( ( t[2] & DIRTYFLAG ) == 0 )
continue;
2523 r = t + t[1]; m = t + FUNHEAD;
2525 for ( j = i+1; j < nden; j++ ) pden[j-1] = pden[j];
2527 for ( j = 0; j < nnco; j++ )
if ( pnco[j] == t )
break;
2528 for ( j++; j < nnco; j++ ) pnco[j-1] = pnco[j];
2534 if ( m >= r )
continue;
2539 k = 1; to = termout; from = term;
2541 while ( from < t ) *to++ = *from++;
2545 *to++ = DENOMINATOR;
2546 for ( j = 1; j < FUNHEAD; j++ ) *to++ = 0;
2547 if ( *m < -FUNCTION ) *to++ = *m++;
2548 else if ( *m < 0 ) { *to++ = *m++; *to++ = *m++; }
2550 j = *m;
while ( --j >= 0 ) *to++ = *m++;
2552 stop[1] = WORDDIF(to,stop);
2555 if ( i == nden - 1 ) {
2556 stop = term + *term;
2557 while ( from < stop ) *to++ = *from++;
2558 i = *termout = WORDDIF(to,termout);
2559 to = term; from = termout;
2560 while ( --i >= 0 ) *to++ = *from++;
2565 for ( i = 0; i < nden; i++ ) {
2567 if ( ( t[2] & DIRTYFLAG ) == 0 )
continue;
2569 if ( t[FUNHEAD] == -SYMBOL ) {
2572 change = ExtraSymbol(*t,-1,nsym,ppsym,&ncoef);
2574 ppsym += change << 1;
2577 else if ( t[FUNHEAD] == -SNUMBER ) {
2579 if ( *t == 0 )
goto NormInf;
2580 if ( *t < 0 ) { *AT.WorkPointer = -*t; j = -1; }
2581 else { *AT.WorkPointer = *t; j = 1; }
2582 ncoef = REDLENG(ncoef);
2583 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,j) )
2585 ncoef = INCLENG(ncoef);
2588 else if ( t[FUNHEAD] == ARGHEAD )
goto NormInf;
2589 else if ( t[FUNHEAD] > 0 && t[FUNHEAD+ARGHEAD] ==
2590 t[FUNHEAD]-ARGHEAD ) {
2593 t += FUNHEAD + ARGHEAD + 1;
2595 m = r - ABS(*r) + 1;
2596 if ( j != 3 || ( ( *m != 1 ) || ( m[1] != 1 ) ) ) {
2597 ncoef = REDLENG(ncoef);
2598 if ( DivRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,REDLENG(j),(UWORD *)n_coef,&ncoef) )
goto FromNorm;
2599 ncoef = INCLENG(ncoef);
2601 t[-FUNHEAD-ARGHEAD] -= j;
2609 if ( *t == SYMBOL || *t == DOTPRODUCT ) {
2612 pden[i][FUNHEAD] -= k;
2613 pden[i][FUNHEAD+ARGHEAD] -= k;
2618 if ( *t == SYMBOL ) {
2622 change = ExtraSymbol(*t,-t[1],nsym,ppsym,&ncoef);
2624 ppsym += change << 1;
2637 while ( to < stop ) *to++ = *from++;
2642 if ( pden[i][1] == 4+FUNHEAD+ARGHEAD ) {
2644 for ( j = 0; j < nnco; j++ ) {
2645 if ( pden[i] == pnco[j] ) {
2647 while ( j < nnco ) {
2648 pnco[j] = pnco[j+1];
2654 pden[i--] = pden[--nden];
2665 for ( i = 0; i < ndel; i += 2 ) {
2666 if ( t[0] == t[1] ) {
2667 if ( t[0] == EMPTYINDEX ) {}
2668 else if ( *t < AM.OffsetIndex ) {
2669 k = AC.FixIndices[*t];
2670 if ( k < 0 ) { j = -1; k = -k; }
2671 else if ( k > 0 ) j = 1;
2675 else if ( *t >= AM.DumInd ) {
2677 if ( k )
goto docontract;
2679 else if ( *t >= AM.WilInd ) {
2680 k = indices[*t-AM.OffsetIndex-WILDOFFSET].dimension;
2681 if ( k )
goto docontract;
2683 else if ( ( k = indices[*t-AM.OffsetIndex].dimension ) != 0 ) {
2687 WithFix: shortnum = k;
2688 ncoef = REDLENG(ncoef);
2689 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),j) )
2691 ncoef = INCLENG(ncoef);
2695 change = ExtraSymbol((WORD)(-k),(WORD)1,nsym,ppsym,&ncoef);
2697 ppsym += change << 1;
2699 t[1] = pdel[ndel-1];
2700 t[0] = pdel[ndel-2];
2707 if ( *t < AM.OffsetIndex && t[1] < AM.OffsetIndex )
goto NormZero;
2708 j = *t - AM.OffsetIndex;
2709 if ( j >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2710 || ( *t < AM.WilInd && indices[j].dimension ) ) ) {
2711 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2714 *m++ = pdel[ndel-2];
2718 else if ( *t == m[1] ) {
2720 *m++ = pdel[ndel-2];
2726 j = t[1]-AM.OffsetIndex;
2727 if ( j >= 0 && ( ( t[1] >= AM.DumInd && AC.lDefDim )
2728 || ( t[1] < AM.WilInd && indices[j].dimension ) ) ) {
2729 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2732 *m++ = pdel[ndel-2];
2736 else if ( t[1] == m[1] ) {
2738 *m++ = pdel[ndel-2];
2750 for ( i = 0; i < ndel; i++ ) {
2751 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2752 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2754 for ( j = 1; j < nvec; j += 2 ) {
2758 *t-- = pdel[--ndel];
2763 t[1] = pdel[--ndel];
2766 *t-- = pdel[--ndel];
2775 if ( ndel > 0 && ncon ) {
2777 for ( i = 0; i < ndel; i++ ) {
2778 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2779 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2780 for ( j = 0; j < ncon; j++ ) {
2781 if ( *pcon[j] == *t ) {
2784 *t-- = pdel[--ndel];
2789 t[1] = pdel[--ndel];
2792 *t-- = pdel[--ndel];
2795 for ( j = 0; j < nnco; j++ ) {
2797 if ( r > m && r < m+m[1] ) {
2798 m[2] |= DIRTYSYMFLAG;
2802 for ( j = 0; j < ncom; j++ ) {
2804 if ( r > m && r < m+m[1] ) {
2805 m[2] |= DIRTYSYMFLAG;
2809 for ( j = 0; j < neps; j++ ) {
2811 if ( r > m && r < m+m[1] ) {
2812 m[2] |= DIRTYSYMFLAG;
2827 for ( i = 3; i < nvec; i += 2 ) {
2828 k = *t - AM.OffsetIndex;
2829 if ( k >= 0 && ( ( *t > AM.DumInd && AC.lDefDim )
2830 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2832 for ( j = i; j < nvec; j += 2 ) {
2838 *r-- = pvec[--nvec];
2840 *t-- = pvec[--nvec];
2841 *t-- = pvec[--nvec];
2850 if ( nvec > 0 && ncon ) {
2852 for ( i = 1; i < nvec; i += 2 ) {
2853 k = *t - AM.OffsetIndex;
2854 if ( k >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2855 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2856 for ( j = 0; j < ncon; j++ ) {
2857 if ( *pcon[j] == *t ) {
2859 *t-- = pvec[--nvec];
2860 *t-- = pvec[--nvec];
2862 pcon[j] = pcon[--ncon];
2864 for ( j = 0; j < nnco; j++ ) {
2866 if ( r > m && r < m+m[1] ) {
2867 m[2] |= DIRTYSYMFLAG;
2871 for ( j = 0; j < ncom; j++ ) {
2873 if ( r > m && r < m+m[1] ) {
2874 m[2] |= DIRTYSYMFLAG;
2878 for ( j = 0; j < neps; j++ ) {
2880 if ( r > m && r < m+m[1] ) {
2881 m[2] |= DIRTYSYMFLAG;
2899 for ( i = 0; i < nnco; i++ ) {
2901 if ( ( *t >= (FUNCTION+WILDOFFSET)
2902 && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
2903 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2904 && functions[*t-FUNCTION].spec == 0 ) ) {
2910 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
2913 pnco[i][2] |= DIRTYSYMFLAG;
2925 for ( i = 0; i < nnco; i++ ) {
2927 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) && *t != DOLLAREXPRESSION ) {
2929 if ( ( *t >= (FUNCTION+WILDOFFSET)
2930 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
2931 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2932 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
2933 if ( *t >= (FUNCTION+WILDOFFSET) ) {
2935 j = FullSymmetrize(BHEAD t,l);
2938 else j = FullSymmetrize(BHEAD t,l);
2939 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
2940 if ( ( j & 2 ) != 0 )
goto NormZero;
2941 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
2944 else t[2] &= ~DIRTYSYMFLAG;
2952 for ( i = 0; i < k; i++ ) {
2954 while ( Commute(pnco[j],pnco[j+1]) ) {
2955 t = pnco[j]; pnco[j] = pnco[j+1]; pnco[j+1] = t;
2957 while ( l >= 0 && Commute(pnco[l],pnco[l+1]) ) {
2958 t = pnco[l]; pnco[l] = pnco[l+1]; pnco[l+1] = t;
2961 if ( ++j >= k )
break;
2968 for ( i = 0; i < nnco; i++ ) {
2970 if ( *t == IDFUNCTION ) AN.idfunctionflag = 1;
2971 if ( *t >= GAMMA && *t <= GAMMASEVEN ) {
2977 *m++ = stype = t[FUNHEAD];
2982 if ( *t == GAMMAFIVE ) {
2983 gtype = GAMMA5; t += FUNHEAD;
goto onegammamatrix; }
2984 else if ( *t == GAMMASIX ) {
2985 gtype = GAMMA6; t += FUNHEAD;
goto onegammamatrix; }
2986 else if ( *t == GAMMASEVEN ) {
2987 gtype = GAMMA7; t += FUNHEAD;
goto onegammamatrix; }
2992 if ( gtype == GAMMA5 ) {
2993 if ( j == GAMMA1 ) j = GAMMA5;
2994 else if ( j == GAMMA5 ) j = GAMMA1;
2995 else if ( j == GAMMA7 ) ncoef = -ncoef;
2996 if ( nnum & 1 ) ncoef = -ncoef;
2998 else if ( gtype == GAMMA6 || gtype == GAMMA7 ) {
3000 if ( gtype == GAMMA6 ) gtype = GAMMA7;
3001 else gtype = GAMMA6;
3003 if ( j == GAMMA1 ) j = gtype;
3004 else if ( j == GAMMA5 ) {
3006 if ( j == GAMMA7 ) ncoef = -ncoef;
3008 else if ( j != gtype )
goto NormZero;
3011 ncoef = REDLENG(ncoef);
3012 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),1) )
goto FromNorm;
3013 ncoef = INCLENG(ncoef);
3017 *m++ = gtype; nnum++;
3022 }
while ( ( ++i < nnco ) && ( *(t = pnco[i]) >= GAMMA
3023 && *t <= GAMMASEVEN ) && ( t[FUNHEAD] == stype ) );
3026 k = WORDDIF(m,to) - FUNHEAD-1;
3029 while ( --k >= 0 ) *from-- = *--r;
3032 to[1] = WORDDIF(m,to);
3034 else if ( *t < 0 ) {
3035 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3039 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3040 && *t != REPLACEMENT && *t != DOLLAREXPRESSION
3041 && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3053 for ( i = 0; i < ncom; i++ ) {
3055 if ( ( *t >= (FUNCTION+WILDOFFSET)
3056 && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
3057 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3058 && functions[*t-FUNCTION].spec == 0 ) ) {
3064 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3067 pcom[i][2] |= DIRTYSYMFLAG;
3077 for ( i = 0; i < ncom; i++ ) {
3079 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) ) {
3081 if ( ( *t >= (FUNCTION+WILDOFFSET)
3082 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3083 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3084 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3085 if ( *t >= (FUNCTION+WILDOFFSET) ) {
3087 j = FullSymmetrize(BHEAD t,l);
3090 else j = FullSymmetrize(BHEAD t,l);
3091 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3092 if ( ( j & 2 ) != 0 )
goto NormZero;
3093 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3096 else t[2] &= ~DIRTYSYMFLAG;
3105 for ( i = 1; i < ncom; i++ ) {
3106 for ( j = i; j > 0; j-- ) {
3112 if ( *r < 0 ) {
if ( *t >= *r )
goto NextI; }
3113 else {
if ( -*t <= *r )
goto NextI; }
3116 else if ( *r < 0 ) {
3117 if ( *t < -*r )
goto NextI;
3120 else if ( *t != *r ) {
3121 if ( *t < *r )
goto NextI;
3122 jexch: t = pcom[j]; pcom[j] = pcom[jj]; pcom[jj] = t;
3125 if ( AC.properorderflag ) {
3126 if ( ( *t >= (FUNCTION+WILDOFFSET)
3127 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
3128 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3129 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) ) {}
3131 WORD *s1, *s2, *ss1, *ss2;
3132 s1 = t+FUNHEAD; s2 = r+FUNHEAD;
3133 ss1 = t + t[1]; ss2 = r + r[1];
3134 while ( s1 < ss1 && s2 < ss2 ) {
3136 if ( k > 0 )
goto jexch;
3137 if ( k < 0 )
goto NextI;
3141 if ( s1 < ss1 )
goto jexch;
3145 kk = r[1] - FUNHEAD;
3148 while ( k > 0 && kk > 0 ) {
3149 if ( *t < *r )
goto NextI;
3150 else if ( *t++ > *r++ )
goto jexch;
3153 if ( k > 0 )
goto jexch;
3159 kk = r[1] - FUNHEAD;
3162 while ( k > 0 && kk > 0 ) {
3163 if ( *t < *r )
goto NextI;
3164 else if ( *t++ > *r++ )
goto jexch;
3167 if ( k > 0 )
goto jexch;
3173 for ( i = 0; i < ncom; i++ ) {
3175 if ( *t == THETA || *t == THETA2 ) {
3176 if ( ( k = DoTheta(BHEAD t) ) == 0 )
goto NormZero;
3182 else if ( *t == DELTA2 || *t == DELTAP ) {
3183 if ( ( k = DoDelta(t) ) == 0 )
goto NormZero;
3189 else if ( *t == AR.PolyFunInv && AR.PolyFunType == 2 ) {
3194 WORD *mm, *tt = t, numt = 0;
3196 while ( tt < t+t[1] ) { numt++; NEXTARG(tt) }
3198 tt = t; mm = m; k = t[1];
3204 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3205 else { *tt++ = *mm++; *tt++ = *mm++; }
3208 k = *mm; NCOPY(tt,mm,k)
3212 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3213 else { *tt++ = *mm++; *tt++ = *mm++; }
3216 k = *mm; NCOPY(tt,mm,k)
3219 t[2] |= MUSTCLEANPRF;
3223 else if ( *t == AR.PolyFun ) {
3224 if ( AR.PolyFunType == 1 ) {
3225 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3226 t[1] == FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER )
goto NormZero;
3227 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) {
3228 if ( AN.PolyNormFlag == 0 ) {
3229 AN.PolyNormFlag = 1;
3236 else if ( AR.PolyFunType == 2 ) {
3247 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3248 t[1] > FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) {
3249 u = t + FUNHEAD + 2;
3251 if ( *u <= -FUNCTION ) {}
3252 else if ( t[1] == FUNHEAD+4 && t[FUNHEAD+2] == -SNUMBER
3253 && t[FUNHEAD+3] == 0 )
goto NormPRF;
3254 else if ( t[1] == FUNHEAD+4 )
goto NormZero;
3256 else if ( t[1] == *u+FUNHEAD+2 )
goto NormZero;
3259 u = t+FUNHEAD; NEXTARG(u);
3260 if ( *u == -SNUMBER && u[1] == 0 )
goto NormInf;
3262 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3263 else if ( i < ncom-1 && pcom[i+1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3265 if ( AN.PolyNormFlag ) {
3266 if ( AR.PolyFunExp == 0 ) {
3270 else if ( AR.PolyFunExp == 1 ) {
3271 if ( PolyFunMode == 0 ) {
3278 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3284 if ( PolyFunMode == 0 ) {
3291 if ( ExpandRat(BHEAD mmm) != 0 )
3298 if ( AR.PolyFunExp == 0 ) {
3302 else if ( AR.PolyFunExp == 1 ) {
3305 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3312 if ( ExpandRat(BHEAD mmm) != 0 )
3319 else if ( *t > 0 ) {
3320 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3321 && *t != REPLACEMENT && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3326 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3335 if ( ReplaceVeto < 0 ) {
3345 WORD *ma = fillsetexp, *mb, *mc;
3348 if ( *ma != REPLACEMENT ) {
3352 if ( *ma == REPLACEMENT && ReplaceType == -1 ) {
3355 if ( AN.RSsize < 2*ma[1]+SUBEXPSIZE ) {
3356 if ( AN.ReplaceScrat ) M_free(AN.ReplaceScrat,
"AN.ReplaceScrat");
3357 AN.RSsize = 2*ma[1]+SUBEXPSIZE+40;
3358 AN.ReplaceScrat = (WORD *)Malloc1((AN.RSsize+1)*
sizeof(WORD),
"AN.ReplaceScrat");
3361 ReplaceSub = AN.ReplaceScrat;
3362 ReplaceSub += SUBEXPSIZE;
3364 if ( *ma > 0 )
goto NoRep;
3365 if ( *ma <= -FUNCTION ) {
3366 *ReplaceSub++ = FUNTOFUN;
3368 *ReplaceSub++ = -*ma++;
3369 if ( *ma > -FUNCTION )
goto NoRep;
3370 *ReplaceSub++ = -*ma++;
3372 else if ( ma+4 > mb )
goto NoRep;
3374 if ( *ma == -SYMBOL ) {
3375 if ( ma[2] == -SYMBOL && ma+4 <= mb )
3376 *ReplaceSub++ = SYMTOSYM;
3377 else if ( ma[2] == -SNUMBER && ma+4 <= mb ) {
3378 *ReplaceSub++ = SYMTONUM;
3379 if ( ReplaceType == 0 ) {
3380 oldtoprhs = C->numrhs;
3385 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3386 *ReplaceSub++ = SYMTONUM;
3388 *ReplaceSub++ = ma[1];
3397 else if ( ma[2] > 0 ) {
3398 WORD *sstop, *ttstop, n;
3402 while ( ss < sstop ) {
3404 ttstop = tt - ABS(tt[-1]);
3406 while ( ss < ttstop ) {
3407 if ( *ss == INDEX )
goto NoRep;
3413 if ( ReplaceType == 0 ) {
3414 oldtoprhs = C->numrhs;
3418 ss =
AddRHS(AT.ebufnum,1);
3423 while ( --n >= 0 ) *ss++ = *tt++;
3425 C->
rhs[C->numrhs+1] = ss;
3427 *ReplaceSub++ = subtype;
3429 *ReplaceSub++ = ma[1];
3430 *ReplaceSub++ = C->numrhs;
3436 else if ( ( *ma == -VECTOR || *ma == -MINVECTOR ) && ma+4 <= mb ) {
3437 if ( ma[2] == -VECTOR ) {
3438 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOVEC;
3439 else *ReplaceSub++ = VECTOMIN;
3441 else if ( ma[2] == -MINVECTOR ) {
3442 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOMIN;
3443 else *ReplaceSub++ = VECTOVEC;
3449 else if ( ma[2] > 0 ) {
3450 WORD *sstop, *ttstop, *w, *mm, n, count;
3452 if ( *ma == -MINVECTOR ) {
3456 while ( ss < sstop ) {
3465 while ( ss < sstop ) {
3467 ttstop = tt - ABS(tt[-1]);
3470 while ( ss < ttstop ) {
3471 if ( *ss == INDEX ) {
3472 n = ss[1] - 2; ss += 2;
3473 while ( --n >= 0 ) {
3474 if ( *ss < MINSPEC ) count++;
3480 if ( count != 1 )
goto NoRep;
3484 if ( ReplaceType == 0 ) {
3485 oldtoprhs = C->numrhs;
3489 mm =
AddRHS(AT.ebufnum,1);
3490 *ReplaceSub++ = subtype;
3492 *ReplaceSub++ = ma[1];
3493 *ReplaceSub++ = C->numrhs;
3497 while ( (mm + n + 10) > C->
Top )
3499 while ( --n >= 0 ) *mm++ = *w++;
3501 C->
rhs[C->numrhs+1] = mm;
3503 mm =
AddRHS(AT.ebufnum,1);
3507 while ( (mm + n + 13) > C->
Top )
3510 while ( w < sstop ) {
3511 tt = w + *w; ttstop = tt - ABS(tt[-1]);
3513 while ( w < ttstop ) {
3514 if ( *w != INDEX ) {
3523 while ( --n >= 0 ) {
3524 if ( *w >= MINSPEC ) *mm++ = *w++;
3529 if ( n <= 2 ) mm -= 2;
3538 while ( w < tt ) *mm++ = *w++;
3539 *ss = WORDDIF(mm,ss);
3542 C->
rhs[C->numrhs+1] = mm;
3544 if ( mm > C->
Top ) {
3545 MLOCK(ErrorMessageLock);
3546 MesPrint(
"Internal error in Normalize with extra compiler buffer");
3547 MUNLOCK(ErrorMessageLock);
3555 else if ( *ma == -INDEX ) {
3556 if ( ( ma[2] == -INDEX || ma[2] == -VECTOR )
3558 *ReplaceSub++ = INDTOIND;
3559 else if ( ma[1] >= AM.OffsetIndex ) {
3560 if ( ma[2] == -SNUMBER && ma+4 <= mb
3561 && ma[3] >= 0 && ma[3] < AM.OffsetIndex )
3562 *ReplaceSub++ = INDTOIND;
3563 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3564 *ReplaceSub++ = INDTOIND;
3566 *ReplaceSub++ = ma[1];
3577 *ReplaceSub++ = ma[1];
3578 *ReplaceSub++ = ma[3];
3583 AN.ReplaceScrat[1] = ReplaceSub-AN.ReplaceScrat;
3588 while ( mb < m ) *mc++ = *mb++;
3592 if ( ReplaceType > 0 ) {
3593 C->numrhs = oldtoprhs;
3597 if ( ++ReplaceVeto >= 0 )
break;
3608 for ( i = 0; i < neps; i++ ) {
3610 if ( ( t[2] & DIRTYSYMFLAG ) != DIRTYSYMFLAG )
continue;
3611 t[2] &= ~DIRTYSYMFLAG;
3612 if ( AR.Eside == LHSIDE || AR.Eside == LHSIDEX ) {
3621 if ( *r != FUNNYWILD ) { r++;
continue; }
3622 k = r[1]; u = r + 2;
3625 if ( *u != FUNNYWILD ) ncoef = -ncoef;
3628 tt[-2] = FUNNYWILD; tt[-1] = k; m -= 2;
3632 for ( r = t + 1; r < m; r++ ) {
3633 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3634 else if ( *r == *t )
goto NormZero;
3639 for ( r = t + 2; r < tt; r += 2 ) {
3640 if ( r[1] < t[1] ) {
3641 k = r[1]; r[1] = t[1]; t[1] = k; ncoef = -ncoef; }
3642 else if ( r[1] == t[1] )
goto NormZero;
3651 for ( r = t + 1; r < m; r++ ) {
3652 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3653 else if ( *r == *t )
goto NormZero;
3662 for ( i = 0; i < (neps-1); i++ ) {
3664 for ( j = i+1; j < neps; j++ ) {
3666 if ( t[1] > r[1] ) {
3667 peps[i] = m = r; peps[j] = r = t; t = m;
3669 else if ( t[1] == r[1] ) {
3675 m = peps[j]; peps[j] = t; peps[i] = t = m;
3678 else if ( *r++ > *m++ )
break;
3679 }
while ( --k > 0 );
3684 for ( i = 0; i < neps; i++ ) {
3696 for ( i = 0; i < ndel; i += 2, r += 2 ) {
3697 if ( r[1] < r[0] ) { k = *r; *r = r[1]; r[1] = k; }
3699 for ( i = 2; i < ndel; i += 2, t += 2 ) {
3701 for ( j = i; j < ndel; j += 2 ) {
3702 if ( *r > *t ) { r += 2; }
3703 else if ( *r < *t ) {
3704 k = *r; *r++ = *t; *t++ = k;
3705 k = *r; *r++ = *t; *t-- = k;
3708 if ( *++r < t[1] ) {
3709 k = *r; *r = t[1]; t[1] = k;
3727 for ( i = 0; i < nind; i++ ) {
3729 for ( j = i+1; j < nind; j++ ) {
3731 k = *r; *r = *t; *t = k;
3749 for ( i = 2; i < nvec; i += 2 ) {
3751 for ( j = i; j < nvec; j += 2 ) {
3753 if ( *++r < t[1] ) {
3754 k = *r; *r = t[1]; t[1] = k;
3758 else if ( *r < *t ) {
3759 k = *r; *r++ = *t; *t++ = k;
3760 k = *r; *r++ = *t; *t-- = k;
3780 while ( --i >= 0 ) {
3781 if ( *t > t[1] ) { j = *t; *t = t[1]; t[1] = j; }
3787 while ( t < (m-3) ) {
3791 if ( *++r == *++t ) {
3793 if ( ( *r < MAXPOWER && t[1] < MAXPOWER )
3794 || ( *r > -MAXPOWER && t[1] > -MAXPOWER ) ) {
3797 if ( *t > MAXPOWER || *t < -MAXPOWER ) {
3798 MLOCK(ErrorMessageLock);
3799 MesPrint(
"Exponent of dotproduct out of range: %d",*t);
3800 MUNLOCK(ErrorMessageLock);
3816 else if ( *r < *++t ) {
3817 k = *r; *r++ = *t; *t = k;
3822 else if ( *r < *t ) {
3823 k = *r; *r++ = *t; *t++ = k;
3824 k = *r; *r++ = *t; *t = k;
3827 else { r += 2; t--; }
3829 else if ( *r < *t ) {
3830 k = *r; *r++ = *t; *t++ = k;
3831 k = *r; *r++ = *t; *t++ = k;
3832 k = *r; *r++ = *t; *t = k;
3841 if ( ( i = ndot ) > 0 ) {
3856 *m++ = ( i = nsym ) + 2;
3859 if ( t[1] < (2*MAXPOWER) ) {
3860 if ( t[1] & 1 ) { *m++ = 0; *m++ = 1; }
3862 if ( *++t & 2 ) ncoef = -ncoef;
3866 else if ( *t <= NumSymbols && *t > -2*MAXPOWER ) {
3867 if ( ( ( ( t[1] > symbols[*t].maxpower ) && ( symbols[*t].maxpower < MAXPOWER ) ) ||
3868 ( ( t[1] < symbols[*t].minpower ) && ( symbols[*t].minpower > -MAXPOWER ) ) ) &&
3869 ( t[1] < 2*MAXPOWER ) && ( t[1] > -2*MAXPOWER ) ) {
3870 if ( i <= 2 || t[2] != *t )
goto NormZero;
3872 if ( AN.ncmod == 1 && ( AC.modmode & ALSOPOWERS ) != 0 ) {
3873 if ( AC.cmod[0] == 1 ) t[1] = 0;
3874 else if ( t[1] >= 0 ) t[1] = 1 + (t[1]-1)%(AC.cmod[0]-1);
3876 t[1] = -1 - (-t[1]-1)%(AC.cmod[0]-1);
3877 if ( t[1] < 0 ) t[1] += (AC.cmod[0]-1);
3880 if ( ( t[1] < (2*MAXPOWER) && t[1] >= MAXPOWER )
3881 || ( t[1] > -(2*MAXPOWER) && t[1] <= -MAXPOWER ) ) {
3882 MLOCK(ErrorMessageLock);
3883 MesPrint(
"Exponent out of range: %d",t[1]);
3884 MUNLOCK(ErrorMessageLock);
3887 if ( AT.TrimPower && AR.PolyFunVar == *t && t[1] > AR.PolyFunPow ) {
3894 else { *r -= 2; t += 2; }
3897 *m++ = *t++; *m++ = *t++;
3899 }
while ( (i-=2) > 0 ); }
3900 if ( *r <= 2 ) m = r-1;
3906 stop = (WORD *)(((UBYTE *)(termout)) + AM.MaxTer);
3908 if ( ( m + i ) > stop ) {
3909 MLOCK(ErrorMessageLock);
3910 MesPrint(
"Term too complex during normalization");
3911 MUNLOCK(ErrorMessageLock);
3914 if ( ReplaceType >= 0 ) {
3921 if ( ReplaceType == 0 ) {
3922 AT.WorkPointer = termout+*termout;
3923 WildFill(BHEAD term,termout,AN.ReplaceScrat);
3924 termout = term + *term;
3927 AT.WorkPointer = r = termout + *termout;
3928 WildFill(BHEAD r,termout,AN.ReplaceScrat);
3935 r += *term; r -= ABS(r[-1]);
3938 if ( *m >= FUNCTION && m[1] > FUNHEAD &&
3939 functions[*m-FUNCTION].spec != TENSORFUNCTION )
3952 TermFree(n_llnum,
"n_llnum");
3953 TermFree(n_coef,
"NormCoef");
3972 if ( termout < term + *term && termout >= term ) AT.WorkPointer = term + *term;
3973 else AT.WorkPointer = termout;
3979 TermFree(n_llnum,
"n_llnum");
3980 TermFree(n_coef,
"NormCoef");
3984 MLOCK(ErrorMessageLock);
3985 MesPrint(
"Division by zero during normalization");
3986 MUNLOCK(ErrorMessageLock);
3990 MLOCK(ErrorMessageLock);
3991 MesPrint(
"0^0 during normalization of term");
3992 MUNLOCK(ErrorMessageLock);
3996 MLOCK(ErrorMessageLock);
3997 MesPrint(
"0/0 in polyratfun during normalization of term");
3998 MUNLOCK(ErrorMessageLock);
4003 AT.WorkPointer = termout;
4004 TermFree(n_llnum,
"n_llnum");
4005 TermFree(n_coef,
"NormCoef");
4009 TermFree(n_llnum,
"n_llnum");
4010 TermFree(n_coef,
"NormCoef");
4014 MLOCK(ErrorMessageLock);
4016 MUNLOCK(ErrorMessageLock);
4017 TermFree(n_llnum,
"n_llnum");
4018 TermFree(n_coef,
"NormCoef");
4031 WORD ExtraSymbol(WORD sym, WORD pow, WORD nsym, WORD *ppsym, WORD *ncoef)
4039 if ( pow > 2*MAXPOWER || pow < -2*MAXPOWER
4040 || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
4041 MLOCK(ErrorMessageLock);
4042 MesPrint(
"Illegal wildcard power combination.");
4043 MUNLOCK(ErrorMessageLock);
4048 if ( ( sym <= NumSymbols && sym > -MAXPOWER )
4049 && ( symbols[sym].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
4050 *m %= symbols[sym].maxpower;
4051 if ( *m < 0 ) *m += symbols[sym].maxpower;
4052 if ( ( symbols[sym].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
4053 if ( ( ( symbols[sym].maxpower & 1 ) == 0 ) &&
4054 ( *m >= symbols[sym].maxpower/2 ) ) {
4055 *m -= symbols[sym].maxpower/2; *ncoef = -*ncoef;
4060 if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
4061 MLOCK(ErrorMessageLock);
4062 MesPrint(
"Power overflow during normalization");
4063 MUNLOCK(ErrorMessageLock);
4069 { *m = m[2]; m++; *m = m[2]; m++; i++; }
4074 else if ( sym < *m ) {
4082 { m--; m[2] = *m; m--; m[2] = *m; i++; }
4093 WORD DoTheta(
PHEAD WORD *t)
4096 WORD k, *r1, *r2, *tstop, type;
4097 WORD ia, *ta, *tb, *stopa, *stopb;
4098 if ( AC.BracketNormalize )
return(-1);
4103 if ( k <= FUNHEAD )
return(1);
4106 if ( r1 == tstop ) {
4110 if ( *t == ARGHEAD ) {
4111 if ( type == THETA )
return(1);
4115 if ( *t == -SNUMBER ) {
4116 if ( t[1] < 0 )
return(0);
4118 if ( type == THETA2 && t[1] == 0 )
return(0);
4125 if ( *t == ABS(k)+1+ARGHEAD ) {
4126 if ( k > 0 )
return(1);
4136 if ( r2 < tstop )
return(-1);
4141 if ( *t == -SNUMBER && *r1 == -SNUMBER ) {
4142 if ( t[1] > r1[1] )
return(0);
4143 else if ( t[1] < r1[1] ) {
4146 else if ( type == THETA )
return(1);
4149 else if ( t[1] == 0 && *t == -SNUMBER ) {
4151 else if ( *t < *r1 )
return(1);
4152 else if ( *t > *r1 )
return(0);
4154 else if ( r1[1] == 0 && *r1 == -SNUMBER ) {
4156 else if ( *t < *r1 )
return(1);
4157 else if ( *t > *r1 )
return(0);
4159 r2 = AT.WorkPointer;
4173 ta += ARGHEAD; tb += ARGHEAD;
4174 while ( ta < stopa ) {
4175 if ( tb >= stopb )
return(0);
4176 if ( ( ia = CompareTerms(BHEAD ta,tb,(WORD)1) ) < 0 )
return(0);
4177 if ( ia > 0 )
return(1);
4181 if ( type == THETA )
return(1);
4190 WORD DoDelta(WORD *t)
4192 WORD k, *r1, *r2, *tstop, isnum, isnum2, type = *t;
4193 if ( AC.BracketNormalize )
return(-1);
4195 if ( k <= FUNHEAD )
goto argzero;
4196 if ( k == FUNHEAD+ARGHEAD && t[FUNHEAD] == ARGHEAD )
goto argzero;
4203 if ( *t == -SNUMBER ) { isnum = 1; k = t[1]; }
4209 if ( k == *t-ARGHEAD-1 ) isnum = 1;
4213 if ( r1 >= tstop ) {
4214 if ( !isnum )
return(-1);
4215 if ( k == 0 )
goto argzero;
4220 if ( r2 < tstop )
return(-1);
4222 if ( *r1 == -SNUMBER ) { isnum2 = 1; }
4228 if ( k == *r1-ARGHEAD-1 ) isnum2 = 1;
4231 if ( isnum != isnum2 )
return(-1);
4233 while ( t < tstop && r1 < r2 ) {
4235 if ( !isnum )
return(-1);
4240 if ( t != tstop || r1 != r2 ) {
4241 if ( !isnum )
return(-1);
4245 if ( type == DELTA2 )
return(1);
4248 if ( type == DELTA2 )
return(0);
4257 void DoRevert(WORD *fun, WORD *tmp)
4259 WORD *t, *r, *m, *to, *tt, *mm, i, j;
4264 if ( *r == -REVERSEFUNCTION ) {
4266 while ( mm < to ) *m++ = *mm++;
4269 fun[2] |= DIRTYSYMFLAG;
4271 else if ( *r <= -FUNCTION ) r++;
4273 if ( *r == -INDEX && r[1] < MINSPEC ) *r = -VECTOR;
4278 if ( ( *r > ARGHEAD )
4279 && ( r[ARGHEAD+1] == REVERSEFUNCTION )
4280 && ( *r == (r[ARGHEAD]+ARGHEAD) )
4281 && ( r[ARGHEAD] == (r[ARGHEAD+2]+4) )
4282 && ( *(r+*r-3) == 1 )
4283 && ( *(r+*r-2) == 1 )
4284 && ( *(r+*r-1) == 3 ) ) {
4296 while ( --j >= 0 ) {
4299 while ( --i >= 0 ) {
4306 else if ( *t <= -FUNCTION ) *m++ = *t++;
4307 else { *m++ = *t++; *m++ = *t++; }
4318 fun[1] = WORDDIF(t,fun);
4320 fun[2] |= DIRTYSYMFLAG;
4343 #define MAXNUMBEROFNONCOMTERMS 2 4345 WORD DetCommu(WORD *terms)
4347 WORD *t, *tnext, *tstop;
4349 if ( *terms == 0 )
return(0);
4350 if ( terms[*terms] == 0 )
return(0);
4354 tstop = tnext - ABS(tnext[-1]);
4356 while ( t < tstop ) {
4357 if ( *t >= FUNCTION ) {
4358 if ( functions[*t-FUNCTION].commute ) {
4360 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4364 else if ( *t == SUBEXPRESSION ) {
4365 if ( cbuf[t[4]].CanCommu[t[2]] ) {
4367 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4371 else if ( *t == EXPRESSION ) {
4373 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4376 else if ( *t == DOLLAREXPRESSION ) {
4383 if ( cbuf[AM.dbufnum].CanCommu[t[2]] ) {
4385 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4404 WORD DoesCommu(WORD *term)
4408 if ( *term == 0 )
return(0);
4409 tstop = term + *term;
4410 tstop = tstop - ABS(tstop[-1]);
4412 while ( term < tstop ) {
4413 if ( ( *term >= FUNCTION ) && ( functions[*term-FUNCTION].commute ) ) {
4415 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4430 WORD *PolyNormPoly (
PHEAD WORD *Poly) {
4433 WORD *buffer = AT.WorkPointer;
4435 if (
NewSort(BHEAD0) ) { Terminate(-1); }
4441 AR.CompareRoutine = (
void *)&
Compare1;
4447 if (
EndSort(BHEAD buffer,1) < 0 ) {
4448 AR.CompareRoutine = (
void *)&
Compare1;
4452 while ( *p ) p += *p;
4453 AR.CompareRoutine = (
void *)&
Compare1;
4454 AT.WorkPointer = p + 1;
4481 WORD *EvaluateGcd(
PHEAD WORD *subterm)
4484 WORD *oldworkpointer = AT.WorkPointer, *work1, *work2, *work3;
4485 WORD *t, *tt, *ttt, *t1, *t2, *t3, *t4, *tstop;
4488 WORD *lnum=n_llnum+1;
4489 WORD *num1, *num2, *num3, *den1, *den2, *den3;
4490 WORD sizenum1, sizenum2, sizenum3, sizeden1, sizeden2, sizeden3;
4491 int i, isnumeric = 0, numarg = 0 ;
4497 tt = subterm + subterm[1]; t = subterm + FUNHEAD;
4501 if ( *t == -SNUMBER ) {
4504 MLOCK(ErrorMessageLock);
4505 MesPrint(
"Trying to take the GCD involving a zero term.");
4506 MUNLOCK(ErrorMessageLock);
4510 t1 = subterm + FUNHEAD;
4511 while ( gcdnum > 1 && t1 < tt ) {
4512 if ( *t1 == -SNUMBER ) {
4514 if ( stor == 0 )
goto gcdzero;
4515 if ( GcdLong(BHEAD (UWORD *)&stor,1,(UWORD *)&gcdnum,1,
4516 (UWORD *)lnum,&nnum) )
goto FromGCD;
4521 else if ( *t1 == -SYMBOL )
goto gcdisone;
4522 else if ( *t1 < 0 )
goto gcdillegal;
4528 ct = *ttt; *ttt = 0;
4530 t1 = PolyNormPoly(BHEAD t1+ARGHEAD);
4539 while ( t3 > t2 && *t3 == 0 ) { t3--; i--; }
4540 if ( GcdLong(BHEAD (UWORD *)t2,(WORD)i,(UWORD *)&gcdnum,1,
4541 (UWORD *)lnum,&nnum) ) {
4546 if ( gcdnum == 1 ) {
4553 AT.WorkPointer = oldworkpointer;
4555 if ( gcdnum == 1 )
goto gcdisone;
4556 oldworkpointer[0] = 4;
4557 oldworkpointer[1] = gcdnum;
4558 oldworkpointer[2] = 1;
4559 oldworkpointer[3] = 3;
4560 oldworkpointer[4] = 0;
4561 AT.WorkPointer = oldworkpointer + 5;
4562 return(oldworkpointer);
4564 else if ( *t == -SYMBOL ) {
4565 t1 = subterm + FUNHEAD;
4568 if ( *t1 == -SNUMBER )
goto gcdisone;
4569 if ( *t1 == -SYMBOL ) {
4570 if ( t1[1] != i )
goto gcdisone;
4573 if ( *t1 < 0 )
goto gcdillegal;
4575 ct = *ttt; *ttt = 0;
4577 t2 = PolyNormPoly(BHEAD t1+ARGHEAD);
4579 else t2 = t1 + ARGHEAD;
4583 tstop = t2 - ABS(t2[-1]);
4584 while ( t3 < tstop ) {
4585 if ( *t3 != SYMBOL ) {
4592 if ( *t4 == i && t4[1] > 0 )
goto nextterminarg;
4602 AT.WorkPointer = oldworkpointer;
4604 oldworkpointer[0] = 8;
4605 oldworkpointer[1] = SYMBOL;
4606 oldworkpointer[2] = 4;
4607 oldworkpointer[3] = t[1];
4608 oldworkpointer[4] = 1;
4609 oldworkpointer[5] = 1;
4610 oldworkpointer[6] = 1;
4611 oldworkpointer[7] = 3;
4612 oldworkpointer[8] = 0;
4613 AT.WorkPointer = oldworkpointer+9;
4614 return(oldworkpointer);
4616 else if ( *t < 0 ) {
4618 MLOCK(ErrorMessageLock);
4619 MesPrint(
"Illegal object in gcd_ function. Object not a number or a symbol.");
4620 MUNLOCK(ErrorMessageLock);
4623 else if ( ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4624 else if ( t[1] != 0 ) {
4625 ttt = t + *t; ct = *ttt; *ttt = 0;
4626 t = PolyNormPoly(BHEAD t+ARGHEAD);
4628 if ( t[*t] == 0 && ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4629 AT.WorkPointer = oldworkpointer;
4644 AT.WorkPointer = oldworkpointer;
4646 t = subterm + FUNHEAD;
4647 for ( i = 1; i < isnumeric; i++ ) {
4651 ttt = t + *t; ct = *ttt; *ttt = 0;
4652 t = PolyNormPoly(BHEAD t+ARGHEAD);
4656 i = (ABS(t[-1])-1)/2;
4659 sizenum1 = sizeden1 = i;
4660 while ( sizenum1 > 1 && num1[sizenum1-1] == 0 ) sizenum1--;
4661 while ( sizeden1 > 1 && den1[sizeden1-1] == 0 ) sizeden1--;
4662 work1 = AT.WorkPointer+1; work2 = work1+sizenum1;
4663 for ( i = 0; i < sizenum1; i++ ) work1[i] = num1[i];
4664 for ( i = 0; i < sizeden1; i++ ) work2[i] = den1[i];
4665 num1 = work1; den1 = work2;
4666 AT.WorkPointer = work2 = work2 + sizeden1;
4667 t = subterm + FUNHEAD;
4669 ttt = t + *t; ct = *ttt; *ttt = 0;
4671 t = PolyNormPoly(BHEAD t+ARGHEAD);
4676 i = (ABS(t[-1])-1)/2;
4679 sizenum2 = sizeden2 = i;
4680 while ( sizenum2 > 1 && num2[sizenum2-1] == 0 ) sizenum2--;
4681 while ( sizeden2 > 1 && den2[sizeden2-1] == 0 ) sizeden2--;
4682 num3 = AT.WorkPointer;
4683 if ( GcdLong(BHEAD (UWORD *)num2,sizenum2,(UWORD *)num1,sizenum1,
4684 (UWORD *)num3,&sizenum3) )
goto FromGCD;
4685 sizenum1 = sizenum3;
4686 for ( i = 0; i < sizenum1; i++ ) num1[i] = num3[i];
4687 den3 = AT.WorkPointer;
4688 if ( GcdLong(BHEAD (UWORD *)den2,sizeden2,(UWORD *)den1,sizeden1,
4689 (UWORD *)den3,&sizeden3) )
goto FromGCD;
4690 sizeden1 = sizeden3;
4691 for ( i = 0; i < sizeden1; i++ ) den1[i] = den3[i];
4692 if ( sizenum1 == 1 && num1[0] == 1 && sizeden1 == 1 && den1[1] == 1 )
4697 AT.WorkPointer = work2;
4699 AT.WorkPointer = oldworkpointer;
4703 if ( sizenum1 > sizeden1 ) {
4704 while ( sizenum1 > sizeden1 ) den1[sizeden1++] = 0;
4706 else if ( sizenum1 < sizeden1 ) {
4707 while ( sizenum1 < sizeden1 ) num1[sizenum1++] = 0;
4712 if ( num1 != t ) { NCOPY(t,num1,sizenum1); }
4714 if ( den1 != t ) { NCOPY(t,den1,sizeden1); }
4719 return(oldworkpointer);
4726 t = subterm + FUNHEAD;
4727 AT.WorkPointer += AM.MaxTer/
sizeof(WORD);
4728 work2 = AT.WorkPointer;
4735 work1 = AT.WorkPointer;
4736 ttt = t + *t; ct = *ttt; *ttt = 0;
4737 t = PolyNormPoly(BHEAD t+ARGHEAD);
4738 if ( *work1 < AT.WorkPointer-work1 ) {
4747 *AT.WorkPointer++ = 0;
4753 if ( work2 != work3 ) {
4754 work1 = PolyGCD2(BHEAD work1,work2);
4756 while ( *work2 ) work2 += *work2;
4760 while ( *work2 ) work2 += *work2;
4761 size = work2 - work1 + 1;
4763 NCOPY(t,work1,size);
4765 return(oldworkpointer);
4768 oldworkpointer[0] = 4;
4769 oldworkpointer[1] = 1;
4770 oldworkpointer[2] = 1;
4771 oldworkpointer[3] = 3;
4772 oldworkpointer[4] = 0;
4773 AT.WorkPointer = oldworkpointer+5;
4774 return(oldworkpointer);
4776 MLOCK(ErrorMessageLock);
4777 MesCall(
"EvaluateGcd");
4778 MUNLOCK(ErrorMessageLock);
4795 int TreatPolyRatFun(
PHEAD WORD *prf)
4797 WORD *t, *tstop, *r, *rstop, *m, *mstop;
4798 WORD exp1 = MAXPOWER, exp2 = MAXPOWER;
4801 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
4802 if ( exp1 > 1 ) exp1 = 1;
4806 if ( exp1 > 0 ) exp1 = 0;
4813 while ( t < tstop ) {
4819 rstop = t - ABS(t[-1]);
4820 while ( r < rstop ) {
4821 if ( *r != SYMBOL ) { r += r[1];
continue; }
4825 while ( m < mstop ) {
4826 if ( *m == AR.PolyFunVar ) {
4827 if ( m[1] < exp1 ) exp1 = m[1];
4833 if ( exp1 > 0 ) exp1 = 0;
4838 if ( exp1 > 0 ) exp1 = 0;
4844 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
4845 if ( exp2 > 1 ) exp2 = 1;
4848 if ( exp2 > 0 ) exp2 = 0;
4854 while ( t < tstop ) {
4860 rstop = t - ABS(t[-1]);
4861 while ( r < rstop ) {
4862 if ( *r != SYMBOL ) { r += r[1];
continue; }
4866 while ( m < mstop ) {
4867 if ( *m == AR.PolyFunVar ) {
4868 if ( m[1] < exp2 ) exp2 = m[1];
4874 if ( exp2 > 0 ) exp2 = 0;
4879 if ( exp2 > 0 ) exp2 = 0;
4892 *t++ = -SNUMBER; *t++ = 1;
4893 *t++ = -SNUMBER; *t++ = 1;
4895 else if ( exp1 > 0 ) {
4897 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
4903 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
4904 *t++ = exp1; *t++ = 1; *t++ = 1; *t++ = 3;
4906 *t++ = -SNUMBER; *t++ = 1;
4909 *t++ = -SNUMBER; *t++ = 1;
4911 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
4917 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
4918 *t++ = -exp1; *t++ = 1; *t++ = 1; *t++ = 3;
4931 void DropCoefficient(
PHEAD WORD *term)
4934 WORD *t = term + *term;
4936 n = t[-1]; na = ABS(n);
4938 if ( n == 3 && t[0] == 1 && t[1] == 1 )
return;
4940 t[0] = 1; t[1] = 1; t[2] = 3;
4949 void DropSymbols(
PHEAD WORD *term)
4952 WORD *tend = term + *term, *t1, *t2, *tstop;
4953 tstop = tend - ABS(tend[-1]);
4955 while ( t1 < tstop ) {
4956 if ( *t1 == SYMBOL ) {
4959 while ( t2 < tend ) *t1++ = *t2++;
4982 WORD buffer[7*NORMSIZE], *t, *b, *bb, *tt, *m, *tstop;
4985 *b++ = SYMBOL; *b++ = 2;
4987 tstop = t - ABS(t[-1]);
4989 while ( t < tstop ) {
4990 if ( *t == SYMBOL && t < tstop ) {
4991 for ( i = 2; i < t[1]; i += 2 ) {
4994 if ( bb[0] == t[i] ) {
4996 if ( bb[1] > MAXPOWER || bb[1] < -MAXPOWER ) {
4997 MLOCK(ErrorMessageLock);
4998 MesPrint(
"Power in SymbolNormalize out of range");
4999 MUNLOCK(ErrorMessageLock);
5005 bb[0] = bb[2]; bb[1] = bb[3]; bb += 2;
5010 else if ( bb[0] > t[i] ) {
5012 while ( m > bb ) { m[1] = m[-1]; m[0] = m[-2]; m -= 2; }
5021 *b++ = t[i]; *b++ = t[i+1];
5027 MLOCK(ErrorMessageLock);
5028 MesPrint(
"Illegal term in SymbolNormalize");
5029 MUNLOCK(ErrorMessageLock);
5034 buffer[1] = b - buffer;
5038 if ( AT.LeaveNegative == 0 ) {
5039 b = buffer; bb = b + b[1]; b += 3;
5042 MLOCK(ErrorMessageLock);
5043 MesPrint(
"Negative power in SymbolNormalize");
5044 MUNLOCK(ErrorMessageLock);
5056 b = buffer; tt = term + 1;
5057 if ( i > 2 ) { NCOPY(tt,b,i) }
5060 if ( i < 0 ) i = -i;
5061 *term -= (tstop-tt);
5075 int TestFunFlag(
PHEAD WORD *tfun)
5077 WORD *t, *tstop, *r, *rstop, *m, *mstop;
5078 if ( functions[*tfun-FUNCTION].spec )
return(0);
5079 tstop = tfun + tfun[1];
5081 while ( t < tstop ) {
5082 if ( *t < 0 ) { NEXTARG(t);
continue; }
5084 if ( t[1] == 0 ) { t = rstop;
continue; }
5086 while ( r < rstop ) {
5087 m = r+1; mstop = r+*r; mstop -= ABS(mstop[-1]);
5088 while ( m < mstop ) {
5089 if ( *m == SUBEXPRESSION || *m == EXPRESSION || *m == DOLLAREXPRESSION )
return(1);
5090 if ( ( *m >= FUNCTION ) && ( ( m[2] & DIRTYFLAG ) == DIRTYFLAG )
5091 && ( *m != REPLACEMENT ) && TestFunFlag(BHEAD m) )
return(1);
5106 #define EXCHN(t1,t2,n) { WORD a,i; for(i=0;i<n;i++){x=t1[i];t1[i]=t2[i];t2[i]=a;} 5107 #define EXCH(x,y) { WORD a = (x); (x) = (y); (y) = a; } 5109 WORD BracketNormalize(
PHEAD WORD *term)
5111 WORD *stop = term+*term-3, *t, *tt, *tstart, *r;
5112 WORD *oldwork = AT.WorkPointer;
5115 termout = AT.WorkPointer = term+*term;
5119 tt = termout+1; t = term+1;
5120 while ( t < stop ) {
5121 if ( *t >= FUNCTION ) { i = t[1]; NCOPY(tt,t,i); }
5124 if ( tt > termout+1 && tt-termout-1 > termout[2] ) {
5125 r = termout+1; ii = tt-r;
5126 for ( i = 0; i < ii-FUNHEAD; i += FUNHEAD ) {
5127 for ( j = i+FUNHEAD; j > 0; j -= FUNHEAD ) {
5128 if ( functions[r[j-FUNHEAD]-FUNCTION].commute
5129 && functions[r[j]-FUNCTION].commute == 0 )
break;
5130 if ( r[j-FUNHEAD] > r[j] ) EXCH(r[j-FUNHEAD],r[j])
5136 tstart = tt; t = term + 1; *tt++ = DELTA; *tt++ = 2;
5137 while ( t < stop ) {
5138 if ( *t == DELTA ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5141 if ( tstart[1] > 2 ) {
5142 for ( r = tstart+2; r < tstart+tstart[1]; r += 2 ) {
5143 if ( r[0] > r[1] ) EXCH(r[0],r[1])
5146 if ( tstart[1] > 4 ) {
5147 r = tstart+2; ii = tstart[1]-2;
5148 for ( i = 0; i < ii-2; i += 2 ) {
5149 for ( j = i+2; j > 0; j -= 2 ) {
5150 if ( r[j-2] > r[j] ) {
5154 else if ( r[j-2] < r[j] )
break;
5156 if ( r[j-1] > r[j+1] ) EXCH(r[j-1],r[j+1])
5161 tt = tstart+tstart[1];
5163 else if ( tstart[1] == 2 ) { tt = tstart; }
5166 tstart = tt; t = term + 1; *tt++ = INDEX; *tt++ = 2;
5167 while ( t < stop ) {
5168 if ( *t == INDEX ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5171 if ( tstart[1] >= 4 ) {
5172 r = tstart+2; ii = tstart[1]-2;
5173 for ( i = 0; i < ii-1; i += 1 ) {
5174 for ( j = i+1; j > 0; j -= 1 ) {
5175 if ( r[j-1] > r[j] ) EXCH(r[j-1],r[j])
5179 tt = tstart+tstart[1];
5181 else if ( tstart[1] == 2 ) { tt = tstart; }
5184 tstart = tt; t = term + 1; *tt++ = DOTPRODUCT; *tt++ = 2;
5185 while ( t < stop ) {
5186 if ( *t == DOTPRODUCT ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5189 if ( tstart[1] > 5 ) {
5190 r = tstart+2; ii = tstart[1]-2;
5191 for ( i = 0; i < ii; i += 3 ) {
5192 if ( r[i] > r[i+1] ) EXCH(r[i],r[i+1])
5194 for ( i = 0; i < ii-3; i += 3 ) {
5195 for ( j = i+3; j > 0; j -= 3 ) {
5196 if ( r[j-3] < r[j] )
break;
5197 if ( r[j-3] > r[j] ) {
5202 if ( r[j-2] > r[j+1] ) EXCH(r[j-2],r[j+1])
5207 tt = tstart+tstart[1];
5209 else if ( tstart[1] == 2 ) { tt = tstart; }
5211 if ( tstart[2] > tstart[3] ) EXCH(tstart[2],tstart[3])
5215 tstart = tt; t = term + 1; *tt++ = SYMBOL; *tt++ = 2;
5216 while ( t < stop ) {
5217 if ( *t == SYMBOL ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5220 if ( tstart[1] > 4 ) {
5221 r = tstart+2; ii = tstart[1]-2;
5222 for ( i = 0; i < ii-2; i += 2 ) {
5223 for ( j = i+2; j > 0; j -= 2 ) {
5224 if ( r[j-2] > r[j] ) EXCH(r[j-2],r[j])
5228 tt = tstart+tstart[1];
5230 else if ( tstart[1] == 2 ) { tt = tstart; }
5233 tstart = tt; t = term + 1; *tt++ = SETSET; *tt++ = 2;
5234 while ( t < stop ) {
5235 if ( *t == SETSET ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5238 if ( tstart[1] > 4 ) {
5239 r = tstart+2; ii = tstart[1]-2;
5240 for ( i = 0; i < ii-2; i += 2 ) {
5241 for ( j = i+2; j > 0; j -= 2 ) {
5242 if ( r[j-2] > r[j] ) {
5249 tt = tstart+tstart[1];
5251 else if ( tstart[1] == 2 ) { tt = tstart; }
5253 *tt++ = 1; *tt++ = 1; *tt++ = 3;
5254 t = term; i = *termout = tt - termout; tt = termout;
5256 AT.WorkPointer = oldwork;
int SymbolNormalize(WORD *term)
WORD StoreTerm(PHEAD WORD *)
WORD Compare1(PHEAD WORD *, WORD *, WORD)
WORD NextPrime(PHEAD WORD)
int CompareSymbols(PHEAD WORD *, WORD *, WORD)
WORD CompCoef(WORD *, WORD *)
LONG EndSort(PHEAD WORD *, int)