48WORD 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);
112WORD 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));
187WORD 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, x;
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[0] == AM.vectorzero )
goto NormZero;
466 if ( t[1] == FUNNYVEC ) {
470 else if ( t[1] < 0 ) {
471 if ( *t == NOINDEX && t[1] == NOINDEX ) t += 2;
473 if ( t[1] == AM.vectorzero )
goto NormZero;
474 *ppdot++ = *t++; *ppdot++ = *t++; *ppdot++ = 1; ndot++;
477 else { *ppvec++ = *t++; *ppvec++ = *t++; nvec += 2; }
483 if ( t[2] == 0 ) t += 3;
484 else if ( ndot > 0 && t[0] == ppdot[-3]
485 && t[1] == ppdot[-2] ) {
488 if ( ppdot[-1] == 0 ) { ppdot -= 3; ndot--; }
490 else if ( t[0] == AM.vectorzero || t[1] == AM.vectorzero ) {
491 if ( t[2] > 0 )
goto NormZero;
495 *ppdot++ = *t++; *ppdot++ = *t++;
496 *ppdot++ = *t++; ndot++;
503 if ( WildFill(BHEAD termout,term,AT.dummysubexp) < 0 )
goto FromNorm;
505 t = termout; m = term;
508 case DOLLAREXPRESSION :
515 if ( AR.Eside != LHSIDE ) {
518 int nummodopt, ptype = -1;
519 if ( AS.MultiThreaded ) {
520 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
521 if ( t[2] == ModOptdollars[nummodopt].number )
break;
523 if ( nummodopt < NumModOptdollars ) {
524 ptype = ModOptdollars[nummodopt].type;
525 if ( ptype == MODLOCAL ) {
526 d = ModOptdollars[nummodopt].dstruct+AT.identity;
529 LOCK(d->pthreadslockread);
534 if ( d->type == DOLZERO ) {
536 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
538 if ( t[3] == 0 )
goto NormZZ;
539 if ( t[3] < 0 )
goto NormInf;
542 else if ( d->type == DOLNUMBER ) {
545 nnum = d->where[nnum-1];
546 if ( nnum < 0 ) { ncoef = -ncoef; nnum = -nnum; }
548 for ( i = 1; i <= nnum; i++ ) lnum[i-1] = d->where[i];
550 if ( nnum == 0 || ( nnum == 1 && lnum[0] == 0 ) ) {
552 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
554 if ( t[3] < 0 )
goto NormInf;
555 else if ( t[3] == 0 )
goto NormZZ;
558 if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) )
goto FromNorm;
559 ncoef = REDLENG(ncoef);
561 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) {
563 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
568 else if ( t[3] > 0 ) {
569 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) {
571 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
576 ncoef = INCLENG(ncoef);
578 else if ( d->type == DOLINDEX ) {
579 if ( d->index == 0 ) {
581 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
585 if ( d->index != NOINDEX ) pind[nind++] = d->index;
587 else if ( d->type == DOLTERMS ) {
588 if ( t[3] >= MAXPOWER || t[3] <= -MAXPOWER ) {
589 if ( d->where[0] == 0 )
goto NormZero;
590 if ( d->where[d->where[0]] != 0 ) {
592 MLOCK(ErrorMessageLock);
593 MesPrint(
"!!!Illegal $ expansion with wildcard power!!!");
594 MUNLOCK(ErrorMessageLock);
601 { WORD *td, *tdstop, dj;
603 tdstop = d->where+d->where[0];
604 if ( tdstop[-1] != 3 || tdstop[-2] != 1
605 || tdstop[-3] != 1 )
goto IllDollarExp;
607 if ( td >= tdstop )
goto IllDollarExp;
608 while ( td < tdstop ) {
609 if ( *td == SYMBOL ) {
610 for ( dj = 2; dj < td[1]; dj += 2 ) {
611 if ( td[dj+1] == 1 ) {
616 else if ( td[dj+1] == -1 ) {
621 else goto IllDollarExp;
624 else if ( *td == DOTPRODUCT ) {
625 for ( dj = 2; dj < td[1]; dj += 3 ) {
626 if ( td[dj+2] == 1 ) {
632 else if ( td[dj+2] == -1 ) {
638 else goto IllDollarExp;
641 else goto IllDollarExp;
648 t[0] = SUBEXPRESSION;
652 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
659 if ( *t == DOLLAREXPRESSION ) {
661 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
665 if ( AS.MultiThreaded ) {
666 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
667 if ( t[2] == ModOptdollars[nummodopt].number )
break;
669 if ( nummodopt < NumModOptdollars ) {
670 ptype = ModOptdollars[nummodopt].type;
671 if ( ptype == MODLOCAL ) {
672 d = ModOptdollars[nummodopt].dstruct+AT.identity;
675 LOCK(d->pthreadslockread);
680 if ( d->type == DOLTERMS ) {
681 t[0] = SUBEXPRESSION;
688 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
694 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
696 MLOCK(ErrorMessageLock);
697 MesPrint(
"!!!This $ variation has not been implemented yet!!!");
698 MUNLOCK(ErrorMessageLock);
702 if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
718 if ( *t == SUMMEDIND ) {
719 if ( t[1] < -NMIN4SHIFT ) {
720 k = -t[1]-NMIN4SHIFT;
721 k = ExtraSymbol(k,1,nsym,ppsym,&ncoef);
725 else if ( t[1] == 0 )
goto NormZero;
735 ncoef = REDLENG(ncoef);
736 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
738 ncoef = INCLENG(ncoef);
742 else if ( *t == NOINDEX && t[1] == NOINDEX ) t += 2;
743 else if ( *t == EMPTYINDEX && t[1] == EMPTYINDEX ) {
744 *ppdel++ = *t++; *ppdel++ = *t++; ndel += 2;
748 *ppdot++ = *t++; *ppdot++ = *t++; *ppdot++ = 1; ndot++;
751 *ppvec++ = *t++; *ppvec++ = *t++; nvec += 2;
756 *ppvec++ = t[1]; *ppvec++ = *t; t+=2; nvec += 2;
758 else { *ppdel++ = *t++; *ppdel++ = *t++; ndel += 2; }
766 if ( t[FUNHEAD] == -SNUMBER && t[1] == FUNHEAD+2
767 && t[FUNHEAD+1] >= 0 ) {
768 if ( Factorial(BHEAD t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
771 ncoef = REDLENG(ncoef);
772 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
goto FromNorm;
773 ncoef = INCLENG(ncoef);
775 else pcom[ncom++] = t;
777 case BERNOULLIFUNCTION :
781 if ( ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] >= 0 )
782 && ( t[1] == FUNHEAD+2 || ( t[1] == FUNHEAD+4 &&
783 t[FUNHEAD+2] == -SNUMBER && ABS(t[FUNHEAD+3]) == 1 ) ) ) {
785 if ( Bernoulli(t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
787 if ( nnum == 0 )
goto NormZero;
788 inum = nnum;
if ( inum < 0 ) inum = -inum;
791 while ( lnum[mnum-1] == 0 ) mnum--;
792 if ( nnum < 0 ) mnum = -mnum;
793 ncoef = REDLENG(ncoef);
794 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,mnum) )
goto FromNorm;
796 while ( lnum[inum+mnum-1] == 0 ) mnum--;
797 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(lnum+inum),mnum) )
goto FromNorm;
798 ncoef = INCLENG(ncoef);
799 if ( t[1] == FUNHEAD+4 && t[FUNHEAD+1] == 1
800 && t[FUNHEAD+3] == -1 ) ncoef = -ncoef;
802 else pcom[ncom++] = t;
814 if ( k == 0 )
goto NormZero;
815 *((UWORD *)lnum) = k;
823 if ( *t == -EXPRESSION ) {
824 k = AS.OldNumFactors[t[1]];
826 else if ( *t == -DOLLAREXPRESSION ) {
827 k = Dollars[t[1]].nfactors;
833 if ( k == 0 )
goto NormZero;
834 *((UWORD *)lnum) = k;
841 if ( t[FUNHEAD] < 0 ) {
842 if ( t[FUNHEAD] <= -FUNCTION && t[1] == FUNHEAD+1 )
break;
843 if ( t[FUNHEAD] > -FUNCTION && t[1] == FUNHEAD+2 ) {
844 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] == 0 )
goto NormZero;
850 if ( t[FUNHEAD] > 0 && t[FUNHEAD] == t[1]-FUNHEAD ) {
852 t += FUNHEAD+ARGHEAD;
857 if ( k == 0 )
goto NormZero;
858 *((UWORD *)lnum) = k;
862 else pcom[ncom++] = t;
866 k = CountFun(AN.cTerm,t);
869 k = CountFun(term,t);
871 if ( k == 0 )
goto NormZero;
872 if ( k > 0 ) { *((UWORD *)lnum) = k; nnum = 1; }
873 else { *((UWORD *)lnum) = -k; nnum = -1; }
877 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+2] == -SNUMBER
878 && t[1] == FUNHEAD+4 && t[FUNHEAD+3] > 1 ) {
880 if ( t[FUNHEAD+1] == 0 )
goto NormZero;
881 if ( t[FUNHEAD+1] < 0 ) { t[FUNHEAD+1] = -t[FUNHEAD+1]; sgn = -1; }
883 if ( MakeRational(t[FUNHEAD+1],t[FUNHEAD+3],x1,x1+1) ) {
884 static int warnflag = 1;
886 MesPrint(
"%w Warning: fraction could not be reconstructed in MakeRational_");
889 x1[0] = t[FUNHEAD+1]; x1[1] = 1;
891 if ( sgn < 0 ) { t[FUNHEAD+1] = -t[FUNHEAD+1]; x1[0] = -x1[0]; }
892 if ( x1[0] < 0 ) { sgn = -1; x1[0] = -x1[0]; }
894 ncoef = REDLENG(ncoef);
895 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)x1,sgn,
896 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
897 ncoef = INCLENG(ncoef);
900 WORD narg = 0, *tt, *ttstop, *arg1 = 0, *arg2 = 0;
903 ttstop = t + t[1]; tt = t+FUNHEAD;
904 while ( tt < ttstop ) {
906 if ( narg == 1 ) arg1 = tt;
910 if ( narg != 2 )
goto defaultcase;
911 if ( *arg2 == -SNUMBER && arg2[1] <= 1 )
goto defaultcase;
912 else if ( *arg2 > 0 && ttstop[-1] < 0 )
goto defaultcase;
913 x1 = NumberMalloc(
"Norm-MakeRational");
914 if ( *arg1 == -SNUMBER ) {
915 if ( arg1[1] == 0 ) {
916 NumberFree(x1,
"Norm-MakeRational");
919 if ( arg1[1] < 0 ) { x1[0] = -arg1[1]; nx1 = -1; }
920 else { x1[0] = arg1[1]; nx1 = 1; }
922 else if ( *arg1 > 0 ) {
924 nx1 = (ABS(arg2[-1])-1)/2;
925 tc = arg1+ARGHEAD+1+nx1;
927 NumberFree(x1,
"Norm-MakeRational");
930 for ( i = 1; i < nx1; i++ )
if ( tc[i] != 0 ) {
931 NumberFree(x1,
"Norm-MakeRational");
935 for ( i = 0; i < nx1; i++ ) x1[i] = tc[i];
936 if ( arg2[-1] < 0 ) nx1 = -nx1;
939 NumberFree(x1,
"Norm-MakeRational");
942 x2 = NumberMalloc(
"Norm-MakeRational");
943 if ( *arg2 == -SNUMBER ) {
944 if ( arg2[1] <= 1 ) {
945 NumberFree(x2,
"Norm-MakeRational");
946 NumberFree(x1,
"Norm-MakeRational");
949 else { x2[0] = arg2[1]; nx2 = 1; }
951 else if ( *arg2 > 0 ) {
953 nx2 = (ttstop[-1]-1)/2;
954 tc = arg2+ARGHEAD+1+nx2;
956 NumberFree(x2,
"Norm-MakeRational");
957 NumberFree(x1,
"Norm-MakeRational");
960 for ( i = 1; i < nx2; i++ )
if ( tc[i] != 0 ) {
961 NumberFree(x2,
"Norm-MakeRational");
962 NumberFree(x1,
"Norm-MakeRational");
966 for ( i = 0; i < nx2; i++ ) x2[i] = tc[i];
969 NumberFree(x2,
"Norm-MakeRational");
970 NumberFree(x1,
"Norm-MakeRational");
973 if ( BigLong(x1,ABS(nx1),x2,nx2) >= 0 ) {
974 UWORD *x3 = NumberMalloc(
"Norm-MakeRational");
975 UWORD *x4 = NumberMalloc(
"Norm-MakeRational");
977 DivLong(x1,nx1,x2,nx2,x3,&nx3,x4,&nx4);
978 for ( i = 0; i < ABS(nx4); i++ ) x1[i] = x4[i];
980 NumberFree(x4,
"Norm-MakeRational");
981 NumberFree(x3,
"Norm-MakeRational");
983 xx = (UWORD *)(TermMalloc(
"Norm-MakeRational"));
984 if ( MakeLongRational(BHEAD x1,nx1,x2,nx2,xx,&nxx) ) {
985 static int warnflag = 1;
987 MesPrint(
"%w Warning: fraction could not be reconstructed in MakeRational_");
990 ncoef = REDLENG(ncoef);
991 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,x1,nx1) )
995 ncoef = REDLENG(ncoef);
996 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,xx,nxx,
997 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
999 ncoef = INCLENG(ncoef);
1000 TermFree(xx,
"Norm-MakeRational");
1001 NumberFree(x2,
"Norm-MakeRational");
1002 NumberFree(x1,
"Norm-MakeRational");
1006 if ( t[1] == FUNHEAD && AN.cTerm ) {
1007 ANsr = r; ANsm = m; ANsc = AN.cTerm;
1011 ncoef = REDLENG(ncoef);
1012 nnum = REDLENG(m[-1]);
1014 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,nnum,
1015 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
1016 ncoef = INCLENG(ncoef);
1021 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1022 if ( GetFirstBracket(termout,t[FUNHEAD+1]) < 0 )
goto FromNorm;
1023 if ( *termout == 0 )
goto NormZero;
1024 if ( *termout > 4 ) {
1026 while ( r < m ) *t++ = *r++;
1028 r2 = termout + *termout; r2 -= ABS(r2[-1]);
1029 while ( r < r1 ) *r2++ = *r++;
1031 while ( r3 < r2 ) *t++ = *r3++;
1033 if ( AT.WorkPointer > term && AT.WorkPointer < t )
1041 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1044 POSITION oldondisk = AS.OldOnFile[t[FUNHEAD+1]];
1045 if ( e->replace == NEWLYDEFINEDEXPRESSION ) {
1046 AS.OldOnFile[t[FUNHEAD+1]] = e->onfile;
1048 if ( *t == FIRSTTERM ) {
1049 if ( GetFirstTerm(termout,t[FUNHEAD+1]) < 0 )
goto FromNorm;
1051 else if ( *t == CONTENTTERM ) {
1052 if ( GetContent(termout,t[FUNHEAD+1]) < 0 )
goto FromNorm;
1054 AS.OldOnFile[t[FUNHEAD+1]] = oldondisk;
1055 if ( *termout == 0 )
goto NormZero;
1059 WORD *r1, *r2, *r3, *r4, *r5, nr1, *rterm;
1060 r2 = termout + *termout; lnum = r2 - ABS(r2[-1]);
1061 nnum = REDLENG(r2[-1]);
1063 r1 = term + *term; r3 = r1 - ABS(r1[-1]);
1064 nr1 = REDLENG(r1[-1]);
1065 if ( Mully(BHEAD (UWORD *)lnum,&nnum,(UWORD *)r3,nr1) )
goto FromNorm;
1066 nnum = INCLENG(nnum); nr1 = ABS(nnum); lnum[nr1-1] = nnum;
1067 rterm = TermMalloc(
"FirstTerm/ContentTerm");
1068 r4 = rterm+1; r5 = term+1;
while ( r5 < t ) *r4++ = *r5++;
1069 r5 = termout+1;
while ( r5 < lnum ) *r4++ = *r5++;
1070 r5 = r;
while ( r5 < r3 ) *r4++ = *r5++;
1071 r5 = lnum; NCOPY(r4,r5,nr1);
1073 nr1 = *rterm; r1 = term; r2 = rterm; NCOPY(r1,r2,nr1);
1074 TermFree(rterm,
"FirstTerm/ContentTerm");
1075 if ( AT.WorkPointer > term && AT.WorkPointer < r1 )
1076 AT.WorkPointer = r1;
1080 else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
1081 DOLLARS d = Dollars + t[FUNHEAD+1], newd = 0;
1084 int nummodopt, dtype = -1;
1085 if ( AS.MultiThreaded && ( AC.mparallelflag == PARALLELFLAG ) ) {
1086 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
1087 if ( t[FUNHEAD+1] == ModOptdollars[nummodopt].number )
break;
1089 if ( nummodopt < NumModOptdollars ) {
1090 dtype = ModOptdollars[nummodopt].type;
1091 if ( dtype == MODLOCAL ) {
1092 d = ModOptdollars[nummodopt].dstruct+AT.identity;
1097 if ( d->where && ( d->type == DOLTERMS || d->type == DOLNUMBER ) ) {
1101 if ( ( newd = DolToTerms(BHEAD t[FUNHEAD+1]) ) == 0 )
1104 if ( newd->where[0] == 0 ) {
1105 M_free(newd,
"Copy of dollar variable");
1108 if ( *t == FIRSTTERM ) {
1109 idol = newd->where[0];
1110 for ( ido = 0; ido < idol; ido++ ) termout[ido] = newd->where[ido];
1112 else if ( *t == CONTENTTERM ) {
1114 tterm = newd->where;
1116 for ( ido = 0; ido < idol; ido++ ) termout[ido] = tterm[ido];
1119 if ( ContentMerge(BHEAD termout,tterm) < 0 )
goto FromNorm;
1124 if ( newd->factors ) M_free(newd->factors,
"Dollar factors");
1125 M_free(newd,
"Copy of dollar variable");
1133 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1134 x = TermsInExpression(t[FUNHEAD+1]);
1135multermnum:
if ( x == 0 )
goto NormZero;
1138 if ( x > (LONG)WORDMASK ) { lnum[0] = x & WORDMASK;
1139 lnum[1] = x >> BITSINWORD; nnum = -2;
1141 else { lnum[0] = x; nnum = -1; }
1143 else if ( x > (LONG)WORDMASK ) {
1144 lnum[0] = x & WORDMASK;
1145 lnum[1] = x >> BITSINWORD;
1148 else { lnum[0] = x; nnum = 1; }
1149 ncoef = REDLENG(ncoef);
1150 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1152 ncoef = INCLENG(ncoef);
1154 else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
1155 x = TermsInDollar(t[FUNHEAD+1]);
1158 else { pcom[ncom++] = t; }
1161 case SIZEOFFUNCTION:
1163 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1164 x = SizeOfExpression(t[FUNHEAD+1]);
1167 else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
1168 x = SizeOfDollar(t[FUNHEAD+1]);
1171 else { pcom[ncom++] = t; }
1175 case PATTERNFUNCTION:
1182 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1183 && t[FUNHEAD+1] >= 0 && t[FUNHEAD+2] == -SNUMBER
1184 && t[FUNHEAD+3] >= 0 && t[FUNHEAD+1] >= t[FUNHEAD+3] ) {
1185 if ( t[FUNHEAD+1] > t[FUNHEAD+3] ) {
1186 if ( GetBinom((UWORD *)lnum,&nnum,
1187 t[FUNHEAD+1],t[FUNHEAD+3]) )
goto FromNorm;
1188 if ( nnum == 0 )
goto NormZero;
1192 else pcom[ncom++] = t;
1198 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1199 if ( ( t[FUNHEAD+1] & 1 ) != 0 ) ncoef = -ncoef;
1201 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] )
1202 && ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) ) {
1203 UWORD *numer1,*denom1;
1204 WORD nsize = abs(t[t[1]-1]), nnsize, isize;
1205 nnsize = (nsize-1)/2;
1206 numer1 = (UWORD *)(t + FUNHEAD+ARGHEAD+1);
1207 denom1 = numer1 + nnsize;
1208 for ( isize = 1; isize < nnsize; isize++ ) {
1209 if ( denom1[isize] )
break;
1211 if ( ( denom1[0] != 1 ) || isize < nnsize ) {
1215 if ( ( numer1[0] & 1 ) != 0 ) ncoef = -ncoef;
1229 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1230 if ( t[FUNHEAD+1] < 0 ) ncoef = -ncoef;
1232 else if ( ( t[1] == FUNHEAD+2 ) && ( t[FUNHEAD] == -SYMBOL )
1233 && ( ( t[FUNHEAD+1] <= NumSymbols && t[FUNHEAD+1] > -MAXPOWER )
1234 && ( symbols[t[FUNHEAD+1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) ) {
1244 *m %= symbols[k].maxpower;
1245 if ( ( symbols[k].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
1246 if ( ( ( symbols[k].maxpower & 1 ) == 0 ) &&
1247 ( *m >= symbols[k].maxpower/2 ) ) {
1248 *m -= symbols[k].maxpower/2; ncoef = -ncoef;
1254 { *m = m[2]; m++; *m = m[2]; m++; i++; }
1260 }
while ( k < *m && --i > 0 );
1263 { m--; m[2] = *m; m--; m[2] = *m; i++; }
1271 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] ) ) {
1272 if ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) {
1273 if ( t[t[1]-1] < 0 ) ncoef = -ncoef;
1278 else if ( ( t[FUNHEAD+ARGHEAD]+FUNHEAD+ARGHEAD == t[1] )
1279 && ( t[FUNHEAD+ARGHEAD+1] == SYMBOL ) ) {
1280 WORD *ts = t + FUNHEAD+ARGHEAD+3;
1281 WORD its = ts[-1]-2;
1283 if ( ( *ts != 0 ) && ( ( *ts > NumSymbols || *ts <= -MAXPOWER )
1284 || ( symbols[*ts].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY ) ) {
1293 if ( t[t[1]-1] < 0 ) ncoef = -ncoef;
1294 ts = t + FUNHEAD+ARGHEAD+3;
1305 if ( ( ts[-1] <= NumSymbols && ts[-1] > -MAXPOWER ) &&
1306 ( symbols[ts[-1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
1307 *m %= symbols[ts[-1]].maxpower;
1308 if ( *m < 0 ) *m += symbols[ts[-1]].maxpower;
1309 if ( ( symbols[ts[-1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
1310 if ( ( ( symbols[ts[-1]].maxpower & 1 ) == 0 ) &&
1311 ( *m >= symbols[ts[-1]].maxpower/2 ) ) {
1312 *m -= symbols[ts[-1]].maxpower/2; ncoef = -ncoef;
1319 { *m = m[2]; m++; *m = m[2]; m++; i++; }
1326 }
while ( *ts < *m && --i > 0 );
1329 { m--; m[2] = *m; m--; m[2] = *m; i++; }
1339signogood: pcom[ncom++] = t;
1342 else pcom[ncom++] = t;
1349 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1351 if ( k < 0 ) k = -k;
1352 if ( k == 0 )
goto NormZero;
1353 *((UWORD *)lnum) = k; nnum = 1;
1357 else if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SYMBOL ) {
1359 if ( ( k > NumSymbols || k <= -MAXPOWER )
1360 || ( symbols[k].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
1363 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] )
1364 && ( t[1] == FUNHEAD+ARGHEAD+t[FUNHEAD+ARGHEAD] ) ) {
1365 if ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) {
1367absnosymbols: ts = t + t[1] -1;
1368 ncoef = REDLENG(ncoef);
1369 nnum = REDLENG(*ts);
1370 if ( nnum < 0 ) nnum = -nnum;
1371 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,
1372 (UWORD *)(ts-ABS(*ts)+1),nnum,
1373 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
1374 ncoef = INCLENG(ncoef);
1379 else if ( t[FUNHEAD+ARGHEAD+1] == SYMBOL ) {
1380 WORD *ts = t+FUNHEAD+ARGHEAD+1;
1381 WORD its = ts[1] - 2;
1385 else if ( ( *ts > NumSymbols || *ts <= -MAXPOWER )
1386 || ( symbols[*ts].complex & VARTYPEROOTOFUNITY )
1387 != VARTYPEROOTOFUNITY )
goto absnogood;
1393absnogood: pcom[ncom++] = t;
1396 else pcom[ncom++] = t;
1404 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1405 && t[FUNHEAD+2] == -SNUMBER && t[FUNHEAD+3] > 1 ) {
1407 tmod = (t[FUNHEAD+1]%t[FUNHEAD+3]);
1408 if ( tmod < 0 ) tmod += t[FUNHEAD+3];
1409 if ( *t == MOD2FUNCTION && tmod > t[FUNHEAD+3]/2 )
1410 tmod -= t[FUNHEAD+3];
1412 *((UWORD *)lnum) = -tmod;
1415 else if ( tmod > 0 ) {
1416 *((UWORD *)lnum) = tmod;
1422 else if ( t[1] > t[FUNHEAD+2] && t[FUNHEAD] > 0
1423 && t[FUNHEAD+t[FUNHEAD]] == -SNUMBER
1424 && t[FUNHEAD+t[FUNHEAD]+1] > 1
1425 && t[1] == FUNHEAD+2+t[FUNHEAD] ) {
1426 WORD *ttt = t+FUNHEAD, iii;
1428 if ( *ttt == ttt[ARGHEAD]+ARGHEAD &&
1429 ttt[ARGHEAD] == ABS(iii)+1 ) {
1431 WORD cmod = ttt[*ttt+1];
1433 if ( *t == MODFUNCTION ) {
1434 if ( TakeModulus((UWORD *)(ttt+ARGHEAD+1)
1435 ,&iii,(UWORD *)(&cmod),ncmod,UNPACK|NOINVERSES) )
1439 if ( TakeModulus((UWORD *)(ttt+ARGHEAD+1)
1440 ,&iii,(UWORD *)(&cmod),ncmod,UNPACK|POSNEG|NOINVERSES) )
1443 if ( *t == MOD2FUNCTION && ttt[ARGHEAD+1] > cmod/2 && iii > 0 ) {
1444 ttt[ARGHEAD+1] -= cmod;
1446 if ( ttt[ARGHEAD+1] < 0 ) {
1447 *((UWORD *)lnum) = -ttt[ARGHEAD+1];
1450 else if ( ttt[ARGHEAD+1] > 0 ) {
1451 *((UWORD *)lnum) = ttt[ARGHEAD+1];
1458 else if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1459 *((UWORD *)lnum) = t[FUNHEAD+1];
1460 if ( *lnum == 0 )
goto NormZero;
1464 else if ( ( ( t[FUNHEAD] < 0 ) && ( t[FUNHEAD] == -SNUMBER )
1465 && ( t[1] >= ( FUNHEAD+6+ARGHEAD ) )
1466 && ( t[FUNHEAD+2] >= 4+ARGHEAD )
1467 && ( t[t[1]-1] == t[FUNHEAD+2+ARGHEAD]-1 ) ) ||
1468 ( ( t[FUNHEAD] > 0 )
1469 && ( t[FUNHEAD]-ARGHEAD-1 == ABS(t[FUNHEAD+t[FUNHEAD]-1]) )
1470 && ( t[FUNHEAD+t[FUNHEAD]]-ARGHEAD-1 == t[t[1]-1] ) ) ) {
1474 WORD *ttt = t + t[1], iii, iii1;
1475 UWORD coefbuf[2], *coef2, ncoef2;
1476 iii = (ttt[-1]-1)/2;
1478 if ( ttt[-1] != 1 ) {
1484 for ( iii1 = 0; iii1 < iii; iii1++ ) {
1485 if ( ttt[iii1] != 0 )
goto exitfromhere;
1496 nnum = -1; lnum[0] = -ttt[1]; lnum[1] = 1;
1499 nnum = 1; lnum[0] = ttt[1]; lnum[1] = 1;
1503 nnum = ABS(ttt[ttt[0]-1] - 1);
1504 for ( iii = 0; iii < nnum; iii++ ) {
1505 lnum[iii] = ttt[ARGHEAD+1+iii];
1508 if ( ttt[ttt[0]-1] < 0 ) nnum = -nnum;
1513 ncoef2 = 3; *coef2 = (UWORD)(ttt[1]);
1517 coef2 = (UWORD *)(ttt+ARGHEAD+1);
1518 ncoef2 = (ttt[ttt[0]-1]-1)/2;
1520 if ( TakeModulus((UWORD *)lnum,&nnum,(UWORD *)coef2,ncoef2,
1521 UNPACK|NOINVERSES|FROMFUNCTION) ) {
1524 if ( *t == MOD2FUNCTION && nnum > 0 ) {
1525 UWORD *coef3 = NumberMalloc(
"Mod2Function"), two = 2;
1527 if ( MulLong((UWORD *)lnum,nnum,&two,1,coef3,&ncoef3) )
1529 if ( BigLong(coef3,ncoef3,(UWORD *)coef2,ncoef2) > 0 ) {
1531 AddLong((UWORD *)lnum,nnum,(UWORD *)coef2,ncoef2
1532 ,(UWORD *)lnum,&nnum);
1535 NumberFree(coef3,
"Mod2Function");
1540 ncoef = REDLENG(ncoef);
1541 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1543 ncoef = INCLENG(ncoef);
1545 else pcom[ncom++] = t;
1549 WORD argcount = 0, *tc, *ts, xc, xs, *tcc;
1550 UWORD *Num1, *Num2, *Num3, *Num4;
1551 WORD size1, size2, size3, size4, space;
1552 tc = t+FUNHEAD; ts = t + t[1];
1553 while ( argcount < 3 && tc < ts ) { NEXTARG(tc); argcount++; }
1554 if ( argcount != 2 )
goto defaultcase;
1555 if ( t[FUNHEAD] == -SNUMBER ) {
1556 if ( t[FUNHEAD+1] <= 1 )
goto defaultcase;
1557 if ( t[FUNHEAD+2] == -SNUMBER ) {
1558 if ( t[FUNHEAD+3] <= 1 )
goto defaultcase;
1559 Num2 = NumberMalloc(
"modinverses");
1560 *Num2 = t[FUNHEAD+3]; size2 = 1;
1563 if ( ts[-1] < 0 )
goto defaultcase;
1564 if ( ts[-1] != t[FUNHEAD+2]-ARGHEAD-1 )
goto defaultcase;
1567 if ( *tcc != 1 )
goto defaultcase;
1568 for ( i = 1; i < xs; i++ ) {
1569 if ( tcc[i] != 0 )
goto defaultcase;
1571 Num2 = NumberMalloc(
"modinverses");
1573 for ( i = 0; i < xs; i++ ) Num2[i] = t[FUNHEAD+ARGHEAD+3+i];
1575 Num1 = NumberMalloc(
"modinverses");
1576 *Num1 = t[FUNHEAD+1]; size1 = 1;
1579 tc = t + FUNHEAD + t[FUNHEAD];
1580 if ( tc[-1] < 0 )
goto defaultcase;
1581 if ( tc[-1] != t[FUNHEAD]-ARGHEAD-1 )
goto defaultcase;
1584 if ( *tcc != 1 )
goto defaultcase;
1585 for ( i = 1; i < xc; i++ ) {
1586 if ( tcc[i] != 0 )
goto defaultcase;
1588 if ( *tc == -SNUMBER ) {
1589 if ( tc[1] <= 1 )
goto defaultcase;
1590 Num2 = NumberMalloc(
"modinverses");
1591 *Num2 = tc[1]; size2 = 1;
1594 if ( ts[-1] < 0 )
goto defaultcase;
1595 if ( ts[-1] != t[FUNHEAD+2]-ARGHEAD-1 )
goto defaultcase;
1598 if ( *tcc != 1 )
goto defaultcase;
1599 for ( i = 1; i < xs; i++ ) {
1600 if ( tcc[i] != 0 )
goto defaultcase;
1602 Num2 = NumberMalloc(
"modinverses");
1604 for ( i = 0; i < xs; i++ ) Num2[i] = tc[ARGHEAD+1+i];
1606 Num1 = NumberMalloc(
"modinverses");
1608 for ( i = 0; i < xc; i++ ) Num1[i] = t[FUNHEAD+ARGHEAD+1+i];
1610 Num3 = NumberMalloc(
"modinverses");
1611 Num4 = NumberMalloc(
"modinverses");
1612 GetLongModInverses(BHEAD Num1,size1,Num2,size2
1613 ,Num3,&size3,Num4,&size4);
1622 if ( ( size3 == 1 || size3 == -1 ) && (*Num3&TOPBITONLY) == 0 ) space += 2;
1623 else space += ARGHEAD + 2*ABS(size3) + 2;
1624 if ( ( size4 == 1 || size4 == -1 ) && (*Num4&TOPBITONLY) == 0 ) space += 2;
1625 else space += ARGHEAD + 2*ABS(size4) + 2;
1626 tt = term + *term; u = tt + space;
1627 while ( tt >= ts ) *--u = *--tt;
1628 m += space; r += space;
1631 if ( ( size3 == 1 || size3 == -1 ) && (*Num3&TOPBITONLY) == 0 ) {
1632 *ts++ = -SNUMBER; *ts = (WORD)(*Num3);
1633 if ( size3 < 0 ) *ts = -*ts;
1637 *ts++ = 2*ABS(size3)+ARGHEAD+2;
1638 *ts++ = 0; FILLARG(ts)
1639 *ts++ = 2*ABS(size3)+1;
1640 for ( i = 0; i < ABS(size3); i++ ) *ts++ = Num3[i];
1642 for ( i = 1; i < ABS(size3); i++ ) *ts++ = 0;
1643 if ( size3 < 0 ) *ts++ = 2*size3-1;
1644 else *ts++ = 2*size3+1;
1646 if ( ( size4 == 1 || size4 == -1 ) && (*Num4&TOPBITONLY) == 0 ) {
1647 *ts++ = -SNUMBER; *ts = *Num4;
1648 if ( size4 < 0 ) *ts = -*ts;
1652 *ts++ = 2*ABS(size4)+ARGHEAD+2;
1653 *ts++ = 0; FILLARG(ts)
1654 *ts++ = 2*ABS(size4)+2;
1655 for ( i = 0; i < ABS(size4); i++ ) *ts++ = Num4[i];
1657 for ( i = 1; i < ABS(size4); i++ ) *ts++ = 0;
1658 if ( size4 < 0 ) *ts++ = 2*size4-1;
1659 else *ts++ = 2*size4+1;
1661 NumberFree(Num4,
"modinverses");
1662 NumberFree(Num3,
"modinverses");
1663 NumberFree(Num1,
"modinverses");
1664 NumberFree(Num2,
"modinverses");
1671#ifdef NEWGCDFUNCTION
1677 WORD *num1, *num2, size1, size2, stor1, stor2, *ttt, ti;
1678 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1679 && t[FUNHEAD+2] == -SNUMBER && t[FUNHEAD+1] != 0
1680 && t[FUNHEAD+3] != 0 ) {
1681 stor1 = t[FUNHEAD+1];
1682 stor2 = t[FUNHEAD+3];
1683 if ( stor1 < 0 ) stor1 = -stor1;
1684 if ( stor2 < 0 ) stor2 = -stor2;
1685 num1 = &stor1; num2 = &stor2;
1689 else if ( t[1] > FUNHEAD+4 ) {
1690 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] != 0
1691 && t[FUNHEAD+2] == t[1]-FUNHEAD-2 &&
1692 ABS(t[t[1]-1]) == t[FUNHEAD+2]-1-ARGHEAD ) {
1694 size2 = ABS(num2[-1]);
1697 size2 = (size2-1)/2;
1699 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1700 if ( ti == 1 && ttt[-1] == 1 ) {
1701 stor1 = t[FUNHEAD+1];
1702 if ( stor1 < 0 ) stor1 = -stor1;
1707 else pcom[ncom++] = t;
1709 else if ( t[FUNHEAD] > 0 &&
1710 t[FUNHEAD]-1-ARGHEAD == ABS(t[t[FUNHEAD]+FUNHEAD-1]) ) {
1711 num1 = t + FUNHEAD + t[FUNHEAD];
1712 size1 = ABS(num1[-1]);
1715 size1 = (size1-1)/2;
1717 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1718 if ( ti == 1 && ttt[-1] == 1 ) {
1719 if ( t[1]-FUNHEAD == t[FUNHEAD]+2 && t[t[1]-2] == -SNUMBER
1720 && t[t[1]-1] != 0 ) {
1722 if ( stor2 < 0 ) stor2 = -stor2;
1727 else if ( t[1]-FUNHEAD == t[FUNHEAD]+t[FUNHEAD+t[FUNHEAD]]
1728 && ABS(t[t[1]-1]) == t[FUNHEAD+t[FUNHEAD]] - ARGHEAD-1 ) {
1730 size2 = ABS(num2[-1]);
1733 size2 = (size2-1)/2;
1735 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1736 if ( ti == 1 && ttt[-1] == 1 ) {
1737gcdcalc:
if ( GcdLong(BHEAD (UWORD *)num1,size1,(UWORD *)num2,size2
1738 ,(UWORD *)lnum,&nnum) )
goto FromNorm;
1741 else pcom[ncom++] = t;
1743 else pcom[ncom++] = t;
1745 else pcom[ncom++] = t;
1747 else pcom[ncom++] = t;
1749 else pcom[ncom++] = t;
1753 WORD *gcd = AT.WorkPointer;
1754 if ( ( gcd = EvaluateGcd(BHEAD t) ) == 0 )
goto FromNorm;
1755 if ( *gcd == 4 && gcd[1] == 1 && gcd[2] == 1 && gcd[4] == 0 ) {
1756 AT.WorkPointer = gcd;
1758 else if ( gcd[*gcd] == 0 ) {
1759 WORD *t1, iii, change, *num, *den, numsize, densize;
1760 if ( gcd[*gcd-1] < *gcd-1 ) {
1762 for ( iii = 2; iii < t1[1]; iii += 2 ) {
1763 change = ExtraSymbol(t1[iii],t1[iii+1],nsym,ppsym,&ncoef);
1765 ppsym += change * 2;
1769 iii = t1[-1]; num = t1-iii; numsize = (iii-1)/2;
1770 den = num + numsize; densize = numsize;
1771 while ( numsize > 1 && num[numsize-1] == 0 ) numsize--;
1772 while ( densize > 1 && den[densize-1] == 0 ) densize--;
1773 if ( numsize > 1 || num[0] != 1 ) {
1774 ncoef = REDLENG(ncoef);
1775 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)num,numsize) )
goto FromNorm;
1776 ncoef = INCLENG(ncoef);
1778 if ( densize > 1 || den[0] != 1 ) {
1779 ncoef = REDLENG(ncoef);
1780 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)den,densize) )
goto FromNorm;
1781 ncoef = INCLENG(ncoef);
1783 AT.WorkPointer = gcd;
1797 LONG size = AT.WorkPointer - gcd;
1799 ss =
AddRHS(AT.ebufnum,1);
1803 C->
rhs[C->numrhs+1] = ss;
1806 t[0] = SUBEXPRESSION;
1813 while ( r < tt ) *t++ = *r++;
1822 MesPrint(
" Unexpected call to EvaluateGCD");
1828 if ( t[1] == FUNHEAD )
break;
1830 WORD *ttt = t + FUNHEAD;
1831 WORD *tttstop = t + t[1];
1833 while ( ttt < tttstop ) {
1835 if ( ttt[ARGHEAD]-1 > ABS(ttt[*ttt-1]) )
goto nospec;
1839 if ( *ttt != -SNUMBER )
goto nospec;
1850 for ( iii = 0; iii < ttt[ARGHEAD]; iii++ )
1851 n_llnum[iii] = ttt[ARGHEAD+iii];
1856 if ( ttt[1] == 0 ) {
1857 n_llnum[0] = n_llnum[1] = n_llnum[2] = n_llnum[3] = 0;
1861 if ( ttt[1] > 0 ) { n_llnum[1] = ttt[1]; n_llnum[3] = 3; }
1862 else { n_llnum[1] = -ttt[1]; n_llnum[3] = -3; }
1870 while ( ttt < tttstop ) {
1872 if ( n_llnum[0] == 0 ) {
1873 if ( ( *t == MINFUNCTION && ttt[*ttt-1] < 0 )
1874 || ( *t == MAXFUNCTION && ttt[*ttt-1] > 0 ) )
1880 if ( ( iii > 0 && *t == MINFUNCTION )
1881 || ( iii < 0 && *t == MAXFUNCTION ) ) {
1882 for ( iii = 0; iii < ttt[0]; iii++ )
1883 n_llnum[iii] = ttt[iii];
1889 if ( n_llnum[0] == 0 ) {
1890 if ( ( *t == MINFUNCTION && ttt[1] < 0 )
1891 || ( *t == MAXFUNCTION && ttt[1] > 0 ) )
1894 else if ( ttt[1] == 0 ) {
1895 if ( ( *t == MINFUNCTION && n_llnum[*n_llnum-1] > 0 )
1896 || ( *t == MAXFUNCTION && n_llnum[*n_llnum-1] < 0 ) ) {
1901 tterm[0] = 4; tterm[2] = 1;
1902 if ( ttt[1] < 0 ) { tterm[1] = -ttt[1]; tterm[3] = -3; }
1903 else { tterm[1] = ttt[1]; tterm[3] = 3; }
1905 if ( ( iii > 0 && *t == MINFUNCTION )
1906 || ( iii < 0 && *t == MAXFUNCTION ) ) {
1907 for ( iii = 0; iii < 4; iii++ )
1908 n_llnum[iii] = tterm[iii];
1914 if ( n_llnum[0] == 0 )
goto NormZero;
1915 ncoef = REDLENG(ncoef);
1916 nnum = REDLENG(n_llnum[*n_llnum-1]);
1917 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)lnum,nnum,
1918 (UWORD *)n_coef,&ncoef) )
goto FromNorm;
1919 ncoef = INCLENG(ncoef);
1922 case INVERSEFACTORIAL:
1923 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] >= 0 ) {
1924 if ( Factorial(BHEAD t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
1926 ncoef = REDLENG(ncoef);
1927 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
goto FromNorm;
1928 ncoef = INCLENG(ncoef);
1931nospec: pcom[ncom++] = t;
1935 if ( ( t[FUNHEAD] == -SYMBOL )
1936 && ( t[FUNHEAD+1] > 0 ) && ( t[1] == FUNHEAD+2 ) ) {
1937 *((UWORD *)lnum) = symbols[t[FUNHEAD+1]].maxpower;
1941 else { pcom[ncom++] = t; }
1944 if ( ( t[FUNHEAD] == -SYMBOL )
1945 && ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+2 ) ) {
1946 *((UWORD *)lnum) = symbols[t[FUNHEAD+1]].minpower;
1950 else { pcom[ncom++] = t; }
1953 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER
1954 && t[FUNHEAD+1] > 0 ) {
1955 UWORD xp = (UWORD)(
NextPrime(BHEAD t[FUNHEAD+1]));
1956 ncoef = REDLENG(ncoef);
1957 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,&xp,1) )
goto FromNorm;
1958 ncoef = INCLENG(ncoef);
1960 else goto defaultcase;
1966 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER
1967 && t[FUNHEAD+1] > 0 ) {
1968 WORD val = Moebius(BHEAD t[FUNHEAD+1]);
1969 if ( val == 0 )
goto NormZero;
1970 if ( val < 0 ) ncoef = -ncoef;
1977 ncoef = REDLENG(ncoef);
1978 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(t+3),t[2]) )
goto FromNorm;
1979 ncoef = INCLENG(ncoef);
1984 if ( t[3] & 1 ) ncoef = -ncoef;
1986 else if ( t[2] == 0 ) {
1987 if ( t[3] < 0 )
goto NormInf;
1992 if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) )
goto FromNorm;
1993 ncoef = REDLENG(ncoef);
1995 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1998 else if ( t[3] > 0 ) {
1999 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2002 ncoef = INCLENG(ncoef);
2009 if ( t[1] == FUNHEAD ) {
2010 MLOCK(ErrorMessageLock);
2011 MesPrint(
"Gamma matrix without spin line encountered.");
2012 MUNLOCK(ErrorMessageLock);
2020 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG ) {
2022 t[2] |= DIRTYSYMFLAG;
2025ScanCont:
while ( t < r ) {
2026 if ( *t >= AM.OffsetIndex &&
2027 ( *t >= AM.DumInd || ( *t < AM.WilInd &&
2028 indices[*t-AM.OffsetIndex].dimension ) ) )
2038 if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) )
2040 if ( *rr == -SNUMBER && rr[1] == 1 )
break;
2041 if ( *rr <= -FUNCTION ) k = *rr;
2043 if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) ) {
2044 if ( k == 0 )
goto NormZZ;
2047 if ( *rr == -SNUMBER && rr[1] > 0 && rr[1] < MAXPOWER && k < 0 ) {
2049 if ( functions[k-FUNCTION].commute ) {
2050 for ( i = 0; i < rr[1]; i++ ) pnco[nnco++] = rr-1;
2053 for ( i = 0; i < rr[1]; i++ ) pcom[ncom++] = rr-1;
2057 if ( k == 0 )
goto NormZero;
2058 if ( t[FUNHEAD] == -SYMBOL && *rr == -SNUMBER && t[1] == FUNHEAD+4 ) {
2059 if ( rr[1] < MAXPOWER ) {
2060 t[FUNHEAD+2] = t[FUNHEAD+1]; t += FUNHEAD+2;
2070 t[2] &= ~DIRTYSYMFLAG;
2076 t[2] &= ~DIRTYSYMFLAG;
2083 if ( *t == 0 || *t == AM.vectorzero )
goto NormZero;
2084 if ( *t > 0 && *t < AM.OffsetIndex ) {
2087 ncoef = REDLENG(ncoef);
2088 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2090 ncoef = INCLENG(ncoef);
2092 else if ( *t == NOINDEX ) t++;
2093 else pind[nind++] = *t++;
2096 case SUBEXPRESSION :
2097 if ( t[3] == 0 )
break;
2109 if ( t[2] == 0 )
goto defaultcase;
2110 if ( t[FUNHEAD] != -SNUMBER || t[FUNHEAD+1] < 0 )
goto defaultcase;
2111 if ( t[FUNHEAD+2] == -SNUMBER ) {
2112 if ( t[FUNHEAD+1] == 0 && t[FUNHEAD+3] == 0 )
goto NormZZ;
2113 if ( t[FUNHEAD+1] == 0 )
break;
2114 if ( t[FUNHEAD+3] < 0 ) {
2115 AT.WorkPointer[0] = -t[FUNHEAD+3];
2119 AT.WorkPointer[0] = t[FUNHEAD+3];
2122 AT.WorkPointer[1] = 1;
2124 else if ( t[FUNHEAD+2] == t[1]-FUNHEAD-2
2125 && t[FUNHEAD+2] == t[FUNHEAD+2+ARGHEAD]+ARGHEAD
2126 && ABS(t[t[1]-1]) == t[FUNHEAD+2+ARGHEAD] - 1 ) {
2128 if ( t[FUNHEAD+1] == 0 )
break;
2129 i = t[t[1]-1]; r1 = t + FUNHEAD+ARGHEAD+3;
2132 r2 = AT.WorkPointer;
2133 while ( --i >= 0 ) *r2++ = *r1++;
2135 else goto defaultcase;
2136 if ( TakeRatRoot((UWORD *)AT.WorkPointer,&nc,t[FUNHEAD+1]) ) {
2140 ncoef = REDLENG(ncoef);
2141 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,nc) )
2143 if ( nc < 0 ) nc = -nc;
2144 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(AT.WorkPointer+nc),nc) )
2146 ncoef = INCLENG(ncoef);
2149 case RANDOMFUNCTION :
2151 WORD nnc, nc, nca, nr;
2162 if ( t[1] == FUNHEAD )
goto defaultcase;
2163 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER &&
2164 t[FUNHEAD+1] > 0 ) {
2165 if ( t[FUNHEAD+1] == 1 )
break;
2167 ((UWORD *)AT.WorkPointer)[0] = wranf(BHEAD0);
2168 ((UWORD *)AT.WorkPointer)[1] = wranf(BHEAD0);
2170 if ( ((UWORD *)AT.WorkPointer)[1] == 0 ) {
2172 if ( ((UWORD *)AT.WorkPointer)[0] == 0 ) {
2176 xx = (UWORD)(t[FUNHEAD+1]);
2178 DivLong((UWORD *)AT.WorkPointer,nr
2180 ,((UWORD *)AT.WorkPointer)+4,&nnc
2181 ,((UWORD *)AT.WorkPointer)+2,&nc);
2182 ((UWORD *)AT.WorkPointer)[4] = 0;
2183 ((UWORD *)AT.WorkPointer)[5] = 0;
2184 ((UWORD *)AT.WorkPointer)[6] = 1;
2185 DivLong((UWORD *)AT.WorkPointer+4,3
2187 ,((UWORD *)AT.WorkPointer)+9,&nnc
2188 ,((UWORD *)AT.WorkPointer)+7,&nca);
2189 AddLong((UWORD *)AT.WorkPointer+4,3
2190 ,((UWORD *)AT.WorkPointer)+7,-nca
2191 ,((UWORD *)AT.WorkPointer)+9,&nnc);
2192 if ( BigLong((UWORD *)AT.WorkPointer,nr
2193 ,((UWORD *)AT.WorkPointer)+9,nnc) >= 0 )
goto redoshort;
2197 AT.WorkPointer[2] = (WORD)xx;
2200 ncoef = REDLENG(ncoef);
2201 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+2,nc) )
2203 ncoef = INCLENG(ncoef);
2205 else if ( t[FUNHEAD] > 0 && t[1] == t[FUNHEAD]+FUNHEAD
2206 && ABS(t[t[1]-1]) == t[FUNHEAD]-1-ARGHEAD && t[t[1]-1] > 0 ) {
2207 WORD nna, nnb, nni, nnb2, nnb2a;
2212 nnt = (UWORD *)(t+t[1]-1-nnb);
2213 if ( *nnt != 1 )
goto defaultcase;
2214 for ( nni = 1; nni < nnb; nni++ ) {
2215 if ( nnt[nni] != 0 )
goto defaultcase;
2217 nnt = (UWORD *)(t + FUNHEAD + ARGHEAD + 1);
2219 for ( nni = 0; nni < nnb2; nni++ ) {
2220 ((UWORD *)AT.WorkPointer)[nni] = wranf(BHEAD0);
2223 while ( nnb2a > 0 && ((UWORD *)AT.WorkPointer)[nnb2a-1] == 0 ) nnb2a--;
2225 DivLong((UWORD *)AT.WorkPointer,nnb2a
2227 ,((UWORD *)AT.WorkPointer)+2*nnb2,&nnc
2228 ,((UWORD *)AT.WorkPointer)+nnb2,&nc);
2229 for ( nni = 0; nni < nnb2; nni++ ) {
2230 ((UWORD *)AT.WorkPointer)[nni+2*nnb2] = 0;
2232 ((UWORD *)AT.WorkPointer)[3*nnb2] = 1;
2233 DivLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2235 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc
2236 ,((UWORD *)AT.WorkPointer)+3*nnb2+1,&nca);
2237 AddLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2238 ,((UWORD *)AT.WorkPointer)+3*nnb2+1,-nca
2239 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc);
2240 if ( BigLong((UWORD *)AT.WorkPointer,nnb2a
2241 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,nnc) >= 0 )
goto redoshort;
2245 for ( nni = 0; nni < nnb; nni++ ) {
2246 ((UWORD *)AT.WorkPointer)[nnb2+nni] = nnt[nni];
2250 ncoef = REDLENG(ncoef);
2251 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+nnb2,nc) )
2253 ncoef = INCLENG(ncoef);
2255 else goto defaultcase;
2259 if ( *t == RANPERM && t[1] > FUNHEAD && t[FUNHEAD] <= -FUNCTION ) {
2261 WORD *mm, *ww, *ow = AT.WorkPointer;
2262 WORD *Array, *targ, *argstop, narg = 0, itot;
2266 while ( targ < argstop ) {
2267 narg++; NEXTARG(targ);
2269 WantAddPointers(narg);
2270 pwork = AT.pWorkSpace + AT.pWorkPointer;
2271 targ = t+FUNHEAD+1; narg = 0;
2272 while ( targ < argstop ) {
2273 pwork[narg++] = targ;
2280 ow = AT.WorkPointer;
2281 Array = AT.WorkPointer;
2282 AT.WorkPointer += narg;
2283 for ( i = 0; i < narg; i++ ) Array[i] = i;
2284 for ( i = 2; i <= narg; i++ ) {
2285 itot = (WORD)(iranf(BHEAD i));
2286 for ( j = 0; j < itot; j++ ) CYCLE1(WORD,Array,i)
2288 mm = AT.WorkPointer;
2289 *mm++ = -t[FUNHEAD];
2291 for ( ie = 2; ie < FUNHEAD; ie++ ) *mm++ = t[ie];
2292 for ( i = 0; i < narg; i++ ) {
2293 ww = pwork[Array[i]];
2296 mm = AT.WorkPointer; t++; ww = t;
2297 i = mm[1]; NCOPY(ww,mm,i)
2298 AT.WorkPointer = ow;
2304 if ( ( t[2] & DIRTYFLAG ) != 0 && t[FUNHEAD] <= -FUNCTION
2305 && t[FUNHEAD+1] == -SNUMBER && t[FUNHEAD+2] > 0 ) {
2306 WORD *rr = t+t[1], *mm = t+FUNHEAD+3, *tt, *tt1, *tt2, num = 0;
2310 while ( mm < rr ) { num++; NEXTARG(mm); }
2311 if ( num < t[FUNHEAD+2] ) { pnco[nnco++] = t;
break; }
2312 *t = -t[FUNHEAD]; mm = t+FUNHEAD+3;
2314 while ( --i > 0 ) { NEXTARG(mm); }
2315 tt = TermMalloc(
"Select_");
2318 for ( i = 0; i < *mm; i++ ) *tt1++ = mm[i];
2320 else if ( *mm <= -FUNCTION ) { *tt1++ = *mm; }
2321 else { *tt1++ = mm[0]; *tt1++ = mm[1]; }
2323 while ( tt2 < mm ) *tt1++ = *tt2++;
2324 i = tt1-tt; tt1 = tt; tt2 = t+FUNHEAD;
2326 TermFree(tt,
"Select_");
2328 while ( mm < rr ) *tt2++ = *mm++;
2331 while ( mm < rr ) *tt2++ = *mm++;
2335 else pnco[nnco++] = t;
2346 if ( t[1] <= FUNHEAD )
break;
2351 if ( *to == ARGHEAD )
goto NormZero;
2355 if ( to[ARGHEAD] != j+1 )
goto NoInteg;
2356 if ( rr >= r ) k = -1;
2357 else if ( *rr == ARGHEAD ) { k = 0; rr += ARGHEAD; }
2358 else if ( *rr == -SNUMBER ) { k = rr[1]; rr += 2; }
2360 if ( rr != r )
goto NoInteg;
2361 if ( k > 1 || k < -1 )
goto NoInteg;
2364 i = ( i < 0 ) ? -j: j;
2365 UnPack((UWORD *)to,i,&den,&num);
2373 if ( AN.NoScrat2 == 0 ) {
2374 AN.NoScrat2 = (UWORD *)Malloc1((AM.MaxTal+2)*
sizeof(UWORD),
"Normalize");
2376 if ( DivLong((UWORD *)to,num,(UWORD *)(to+j),den
2377 ,(UWORD *)AT.WorkPointer,&num,AN.NoScrat2,&den) )
goto FromNorm;
2378 if ( k < 0 && den < 0 ) {
2381 if ( AddLong((UWORD *)AT.WorkPointer,num
2382 ,AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2385 else if ( k > 0 && den > 0 ) {
2388 if ( AddLong((UWORD *)AT.WorkPointer,num,
2389 AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2394 else if ( *to == -SNUMBER ) {
2395 if ( to[1] < 0 ) { *AT.WorkPointer = -to[1]; num = -1; }
2396 else if ( to[1] == 0 )
goto NormZero;
2397 else { *AT.WorkPointer = to[1]; num = 1; }
2400 if ( num == 0 )
goto NormZero;
2401 ncoef = REDLENG(ncoef);
2402 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,num) )
2404 ncoef = INCLENG(ncoef);
2413 if ( *t < FUNCTION ) {
2414 MLOCK(ErrorMessageLock);
2415 MesPrint(
"Illegal code in Norm");
2419 AO.OutFill = AO.OutputLine = OutBuf;
2424 while ( --i >= 0 ) {
2425 TalToLine((UWORD)(*t++));
2426 TokenToLine((UBYTE *)
" ");
2432 MUNLOCK(ErrorMessageLock);
2435 if ( *t == REPLACEMENT ) {
2436 if ( AR.Eside != LHSIDE ) ReplaceVeto--;
2444 if ( *t == DUMMYFUN || *t == DUMMYTEN ) {}
2446 if ( *t < (FUNCTION + WILDOFFSET) ) {
2447 if ( ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2448 || ( functions[*t-FUNCTION].minnumargs > 0 ) )
2449 && ( ( t[2] & DIRTYFLAG ) != 0 ) ) {
2453 WORD *ta = t + FUNHEAD, *tb = t + t[1];
2455 while ( ta < tb ) { numarg++; NEXTARG(ta) }
2456 if ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2457 && ( numarg >= functions[*t-FUNCTION].maxnumargs ) )
2459 if ( ( functions[*t-FUNCTION].minnumargs > 0 )
2460 && ( numarg < functions[*t-FUNCTION].minnumargs ) )
2464 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION].tabl == 0 ) ) {
2466 t[2] |= DIRTYSYMFLAG;
2468 if ( functions[*t-FUNCTION].commute ) { pnco[nnco++] = t; }
2469 else { pcom[ncom++] = t; }
2472 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION-WILDOFFSET].tabl == 0 ) ) {
2474 t[2] |= DIRTYSYMFLAG;
2476 if ( functions[*t-FUNCTION-WILDOFFSET].commute ) {
2479 else { pcom[ncom++] = t; }
2485 if ( ( *t < (FUNCTION + WILDOFFSET)
2486 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) || (
2487 *t >= (FUNCTION + WILDOFFSET)
2488 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION ) ) {
2489 if ( *t >= GAMMA && *t <= GAMMASEVEN ) t++;
2492 if ( *t == AM.vectorzero )
goto NormZero;
2493 if ( *t >= AM.OffsetIndex && ( *t >= AM.DumInd
2494 || ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2497 else if ( *t == FUNNYWILD ) { t++; }
2512 else if ( *t <= -FUNCTION ) t++;
2513 else if ( *t == -INDEX ) {
2514 if ( t[1] >= AM.OffsetIndex &&
2515 ( t[1] >= AM.DumInd || ( t[1] < AM.WilInd
2516 && indices[t[1]-AM.OffsetIndex].dimension ) ) )
2520 else if ( *t == -SYMBOL ) {
2521 if ( t[1] >= MAXPOWER && t[1] < 2*MAXPOWER ) {
2525 else if ( t[1] < -MAXPOWER && t[1] > -2*MAXPOWER ) {
2541 r = t = ANsr; m = ANsm;
2542 ANsc = ANsm = ANsr = 0;
2555 for ( k = 0, i = 0; i < nden; i++ ) {
2557 if ( ( t[2] & DIRTYFLAG ) == 0 )
continue;
2558 r = t + t[1]; m = t + FUNHEAD;
2560 for ( j = i+1; j < nden; j++ ) pden[j-1] = pden[j];
2562 for ( j = 0; j < nnco; j++ )
if ( pnco[j] == t )
break;
2563 for ( j++; j < nnco; j++ ) pnco[j-1] = pnco[j];
2569 if ( m >= r )
continue;
2574 k = 1; to = termout; from = term;
2576 while ( from < t ) *to++ = *from++;
2580 *to++ = DENOMINATOR;
2581 for ( j = 1; j < FUNHEAD; j++ ) *to++ = 0;
2582 if ( *m < -FUNCTION ) *to++ = *m++;
2583 else if ( *m < 0 ) { *to++ = *m++; *to++ = *m++; }
2585 j = *m;
while ( --j >= 0 ) *to++ = *m++;
2587 stop[1] = WORDDIF(to,stop);
2590 if ( i == nden - 1 ) {
2591 stop = term + *term;
2592 while ( from < stop ) *to++ = *from++;
2593 i = *termout = WORDDIF(to,termout);
2594 to = term; from = termout;
2595 while ( --i >= 0 ) *to++ = *from++;
2600 for ( i = 0; i < nden; i++ ) {
2602 if ( ( t[2] & DIRTYFLAG ) == 0 )
continue;
2604 if ( t[FUNHEAD] == -SYMBOL ) {
2607 change = ExtraSymbol(*t,-1,nsym,ppsym,&ncoef);
2609 ppsym += change * 2;
2612 else if ( t[FUNHEAD] == -SNUMBER ) {
2614 if ( *t == 0 )
goto NormInf;
2615 if ( *t < 0 ) { *AT.WorkPointer = -*t; j = -1; }
2616 else { *AT.WorkPointer = *t; j = 1; }
2617 ncoef = REDLENG(ncoef);
2618 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,j) )
2620 ncoef = INCLENG(ncoef);
2623 else if ( t[FUNHEAD] == ARGHEAD )
goto NormInf;
2624 else if ( t[FUNHEAD] > 0 && t[FUNHEAD+ARGHEAD] ==
2625 t[FUNHEAD]-ARGHEAD ) {
2628 t += FUNHEAD + ARGHEAD + 1;
2630 m = r - ABS(*r) + 1;
2631 if ( j != 3 || ( ( *m != 1 ) || ( m[1] != 1 ) ) ) {
2632 ncoef = REDLENG(ncoef);
2633 if ( DivRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,REDLENG(j),(UWORD *)n_coef,&ncoef) )
goto FromNorm;
2634 ncoef = INCLENG(ncoef);
2636 t[-FUNHEAD-ARGHEAD] -= j;
2644 if ( *t == SYMBOL || *t == DOTPRODUCT ) {
2647 pden[i][FUNHEAD] -= k;
2648 pden[i][FUNHEAD+ARGHEAD] -= k;
2653 if ( *t == SYMBOL ) {
2657 change = ExtraSymbol(*t,-t[1],nsym,ppsym,&ncoef);
2659 ppsym += change * 2;
2672 while ( to < stop ) *to++ = *from++;
2677 if ( pden[i][1] == 4+FUNHEAD+ARGHEAD ) {
2679 for ( j = 0; j < nnco; j++ ) {
2680 if ( pden[i] == pnco[j] ) {
2682 while ( j < nnco ) {
2683 pnco[j] = pnco[j+1];
2689 pden[i--] = pden[--nden];
2700 for ( i = 0; i < ndel; i += 2 ) {
2701 if ( t[0] == t[1] ) {
2702 if ( t[0] == EMPTYINDEX ) {}
2703 else if ( *t < AM.OffsetIndex ) {
2704 k = AC.FixIndices[*t];
2705 if ( k < 0 ) { j = -1; k = -k; }
2706 else if ( k > 0 ) j = 1;
2710 else if ( *t >= AM.DumInd ) {
2712 if ( k )
goto docontract;
2714 else if ( *t >= AM.WilInd ) {
2715 k = indices[*t-AM.OffsetIndex-WILDOFFSET].dimension;
2716 if ( k )
goto docontract;
2718 else if ( ( k = indices[*t-AM.OffsetIndex].dimension ) != 0 ) {
2722WithFix: shortnum = k;
2723 ncoef = REDLENG(ncoef);
2724 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),j) )
2726 ncoef = INCLENG(ncoef);
2730 change = ExtraSymbol((WORD)(-k),(WORD)1,nsym,ppsym,&ncoef);
2732 ppsym += change * 2;
2734 t[1] = pdel[ndel-1];
2735 t[0] = pdel[ndel-2];
2742 if ( *t < AM.OffsetIndex && t[1] < AM.OffsetIndex )
goto NormZero;
2743 j = *t - AM.OffsetIndex;
2744 if ( j >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2745 || ( *t < AM.WilInd && indices[j].dimension ) ) ) {
2746 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2749 *m++ = pdel[ndel-2];
2753 else if ( *t == m[1] ) {
2755 *m++ = pdel[ndel-2];
2761 j = t[1]-AM.OffsetIndex;
2762 if ( j >= 0 && ( ( t[1] >= AM.DumInd && AC.lDefDim )
2763 || ( t[1] < AM.WilInd && indices[j].dimension ) ) ) {
2764 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2767 *m++ = pdel[ndel-2];
2771 else if ( t[1] == m[1] ) {
2773 *m++ = pdel[ndel-2];
2785 for ( i = 0; i < ndel; i++ ) {
2786 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2787 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2789 for ( j = 1; j < nvec; j += 2 ) {
2793 *t-- = pdel[--ndel];
2798 t[1] = pdel[--ndel];
2801 *t-- = pdel[--ndel];
2810 if ( ndel > 0 && ncon ) {
2812 for ( i = 0; i < ndel; i++ ) {
2813 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2814 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2815 for ( j = 0; j < ncon; j++ ) {
2816 if ( *pcon[j] == *t ) {
2819 *t-- = pdel[--ndel];
2824 t[1] = pdel[--ndel];
2827 *t-- = pdel[--ndel];
2830 for ( j = 0; j < nnco; j++ ) {
2832 if ( r > m && r < m+m[1] ) {
2833 m[2] |= DIRTYSYMFLAG;
2837 for ( j = 0; j < ncom; j++ ) {
2839 if ( r > m && r < m+m[1] ) {
2840 m[2] |= DIRTYSYMFLAG;
2844 for ( j = 0; j < neps; j++ ) {
2846 if ( r > m && r < m+m[1] ) {
2847 m[2] |= DIRTYSYMFLAG;
2862 for ( i = 3; i < nvec; i += 2 ) {
2863 k = *t - AM.OffsetIndex;
2864 if ( k >= 0 && ( ( *t > AM.DumInd && AC.lDefDim )
2865 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2867 for ( j = i; j < nvec; j += 2 ) {
2873 *r-- = pvec[--nvec];
2875 *t-- = pvec[--nvec];
2876 *t-- = pvec[--nvec];
2885 if ( nvec > 0 && ncon ) {
2887 for ( i = 1; i < nvec; i += 2 ) {
2888 k = *t - AM.OffsetIndex;
2889 if ( k >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2890 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2891 for ( j = 0; j < ncon; j++ ) {
2892 if ( *pcon[j] == *t ) {
2894 *t-- = pvec[--nvec];
2895 *t-- = pvec[--nvec];
2897 pcon[j] = pcon[--ncon];
2899 for ( j = 0; j < nnco; j++ ) {
2901 if ( r > m && r < m+m[1] ) {
2902 m[2] |= DIRTYSYMFLAG;
2906 for ( j = 0; j < ncom; j++ ) {
2908 if ( r > m && r < m+m[1] ) {
2909 m[2] |= DIRTYSYMFLAG;
2913 for ( j = 0; j < neps; j++ ) {
2915 if ( r > m && r < m+m[1] ) {
2916 m[2] |= DIRTYSYMFLAG;
2934 for ( i = 0; i < nnco; i++ ) {
2936 if ( ( *t >= (FUNCTION+WILDOFFSET)
2937 && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
2938 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2939 && functions[*t-FUNCTION].spec == 0 ) ) {
2945 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
2948 pnco[i][2] |= DIRTYSYMFLAG;
2960 for ( i = 0; i < nnco; i++ ) {
2962 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) && *t != DOLLAREXPRESSION ) {
2964 if ( ( *t >= (FUNCTION+WILDOFFSET)
2965 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
2966 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2967 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
2968 if ( *t >= (FUNCTION+WILDOFFSET) ) {
2970 j = FullSymmetrize(BHEAD t,l);
2973 else j = FullSymmetrize(BHEAD t,l);
2974 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
2975 if ( ( j & 2 ) != 0 )
goto NormZero;
2976 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
2979 else t[2] &= ~DIRTYSYMFLAG;
2987 for ( i = 0; i < k; i++ ) {
2989 while ( Commute(pnco[j],pnco[j+1]) ) {
2990 t = pnco[j]; pnco[j] = pnco[j+1]; pnco[j+1] = t;
2992 while ( l >= 0 && Commute(pnco[l],pnco[l+1]) ) {
2993 t = pnco[l]; pnco[l] = pnco[l+1]; pnco[l+1] = t;
2996 if ( ++j >= k )
break;
3003 for ( i = 0; i < nnco; i++ ) {
3005 if ( *t == IDFUNCTION ) AN.idfunctionflag = 1;
3006 if ( *t >= GAMMA && *t <= GAMMASEVEN ) {
3012 *m++ = stype = t[FUNHEAD];
3017 if ( *t == GAMMAFIVE ) {
3018 gtype = GAMMA5; t += FUNHEAD;
goto onegammamatrix; }
3019 else if ( *t == GAMMASIX ) {
3020 gtype = GAMMA6; t += FUNHEAD;
goto onegammamatrix; }
3021 else if ( *t == GAMMASEVEN ) {
3022 gtype = GAMMA7; t += FUNHEAD;
goto onegammamatrix; }
3027 if ( gtype == GAMMA5 ) {
3028 if ( j == GAMMA1 ) j = GAMMA5;
3029 else if ( j == GAMMA5 ) j = GAMMA1;
3030 else if ( j == GAMMA7 ) ncoef = -ncoef;
3031 if ( nnum & 1 ) ncoef = -ncoef;
3033 else if ( gtype == GAMMA6 || gtype == GAMMA7 ) {
3035 if ( gtype == GAMMA6 ) gtype = GAMMA7;
3036 else gtype = GAMMA6;
3038 if ( j == GAMMA1 ) j = gtype;
3039 else if ( j == GAMMA5 ) {
3041 if ( j == GAMMA7 ) ncoef = -ncoef;
3043 else if ( j != gtype )
goto NormZero;
3046 ncoef = REDLENG(ncoef);
3047 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),1) )
goto FromNorm;
3048 ncoef = INCLENG(ncoef);
3052 *m++ = gtype; nnum++;
3057 }
while ( ( ++i < nnco ) && ( *(t = pnco[i]) >= GAMMA
3058 && *t <= GAMMASEVEN ) && ( t[FUNHEAD] == stype ) );
3061 k = WORDDIF(m,to) - FUNHEAD-1;
3064 while ( --k >= 0 ) *from-- = *--r;
3067 to[1] = WORDDIF(m,to);
3069 else if ( *t < 0 ) {
3070 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3074 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3075 && *t != REPLACEMENT && *t != DOLLAREXPRESSION
3076 && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3088 for ( i = 0; i < ncom; i++ ) {
3090 if ( ( *t >= (FUNCTION+WILDOFFSET)
3091 && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
3092 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3093 && functions[*t-FUNCTION].spec == 0 ) ) {
3099 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3102 pcom[i][2] |= DIRTYSYMFLAG;
3112 for ( i = 0; i < ncom; i++ ) {
3114 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) ) {
3116 if ( ( *t >= (FUNCTION+WILDOFFSET)
3117 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3118 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3119 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3120 if ( *t >= (FUNCTION+WILDOFFSET) ) {
3122 j = FullSymmetrize(BHEAD t,l);
3125 else j = FullSymmetrize(BHEAD t,l);
3126 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3127 if ( ( j & 2 ) != 0 )
goto NormZero;
3128 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3131 else t[2] &= ~DIRTYSYMFLAG;
3140 for ( i = 1; i < ncom; i++ ) {
3141 for ( j = i; j > 0; j-- ) {
3147 if ( *r < 0 ) {
if ( *t >= *r )
goto NextI; }
3148 else {
if ( -*t <= *r )
goto NextI; }
3151 else if ( *r < 0 ) {
3152 if ( *t < -*r )
goto NextI;
3155 else if ( *t != *r ) {
3156 if ( *t < *r )
goto NextI;
3157jexch: t = pcom[j]; pcom[j] = pcom[jj]; pcom[jj] = t;
3160 if ( AC.properorderflag ) {
3161 if ( ( *t >= (FUNCTION+WILDOFFSET)
3162 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
3163 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3164 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) ) {}
3166 WORD *s1, *s2, *ss1, *ss2;
3167 s1 = t+FUNHEAD; s2 = r+FUNHEAD;
3168 ss1 = t + t[1]; ss2 = r + r[1];
3169 while ( s1 < ss1 && s2 < ss2 ) {
3171 if ( k > 0 )
goto jexch;
3172 if ( k < 0 )
goto NextI;
3176 if ( s1 < ss1 )
goto jexch;
3180 kk = r[1] - FUNHEAD;
3183 while ( k > 0 && kk > 0 ) {
3184 if ( *t < *r )
goto NextI;
3185 else if ( *t++ > *r++ )
goto jexch;
3188 if ( k > 0 )
goto jexch;
3194 kk = r[1] - FUNHEAD;
3197 while ( k > 0 && kk > 0 ) {
3198 if ( *t < *r )
goto NextI;
3199 else if ( *t++ > *r++ )
goto jexch;
3202 if ( k > 0 )
goto jexch;
3208 for ( i = 0; i < ncom; i++ ) {
3210 if ( *t == THETA || *t == THETA2 ) {
3211 if ( ( k = DoTheta(BHEAD t) ) == 0 )
goto NormZero;
3217 else if ( *t == DELTA2 || *t == DELTAP ) {
3218 if ( ( k = DoDelta(t) ) == 0 )
goto NormZero;
3224 else if ( *t == AR.PolyFunInv && AR.PolyFunType == 2 ) {
3229 WORD *mm, *tt = t, numt = 0;
3231 while ( tt < t+t[1] ) { numt++; NEXTARG(tt) }
3233 tt = t; mm = m; k = t[1];
3239 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3240 else { *tt++ = *mm++; *tt++ = *mm++; }
3243 k = *mm; NCOPY(tt,mm,k)
3247 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3248 else { *tt++ = *mm++; *tt++ = *mm++; }
3251 k = *mm; NCOPY(tt,mm,k)
3254 t[2] |= MUSTCLEANPRF;
3258 else if ( *t == AR.PolyFun ) {
3259 if ( AR.PolyFunType == 1 ) {
3260 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3261 t[1] == FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER )
goto NormZero;
3262 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) {
3263 if ( AN.PolyNormFlag == 0 ) {
3264 AN.PolyNormFlag = 1;
3271 else if ( AR.PolyFunType == 2 ) {
3282 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3283 t[1] > FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) {
3284 u = t + FUNHEAD + 2;
3286 if ( *u <= -FUNCTION ) {}
3287 else if ( t[1] == FUNHEAD+4 && t[FUNHEAD+2] == -SNUMBER
3288 && t[FUNHEAD+3] == 0 )
goto NormPRF;
3289 else if ( t[1] == FUNHEAD+4 )
goto NormZero;
3291 else if ( t[1] == *u+FUNHEAD+2 )
goto NormZero;
3294 u = t+FUNHEAD; NEXTARG(u);
3295 if ( *u == -SNUMBER && u[1] == 0 )
goto NormInf;
3297 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3298 else if ( i < ncom-1 && pcom[i+1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3300 if ( AN.PolyNormFlag ) {
3301 if ( AR.PolyFunExp == 0 ) {
3305 else if ( AR.PolyFunExp == 1 ) {
3306 if ( PolyFunMode == 0 ) {
3313 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3319 if ( PolyFunMode == 0 ) {
3326 if ( ExpandRat(BHEAD mmm) != 0 )
3333 if ( AR.PolyFunExp == 0 ) {
3337 else if ( AR.PolyFunExp == 1 ) {
3340 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3347 if ( ExpandRat(BHEAD mmm) != 0 )
3354 else if ( *t > 0 ) {
3355 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3356 && *t != REPLACEMENT && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3361 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3370 if ( ReplaceVeto < 0 ) {
3380 WORD *ma = fillsetexp, *mb, *mc;
3383 if ( *ma != REPLACEMENT ) {
3387 if ( *ma == REPLACEMENT && ReplaceType == -1 ) {
3390 if ( AN.RSsize < 2*ma[1]+SUBEXPSIZE ) {
3391 if ( AN.ReplaceScrat ) M_free(AN.ReplaceScrat,
"AN.ReplaceScrat");
3392 AN.RSsize = 2*ma[1]+SUBEXPSIZE+40;
3393 AN.ReplaceScrat = (WORD *)Malloc1((AN.RSsize+1)*
sizeof(WORD),
"AN.ReplaceScrat");
3396 ReplaceSub = AN.ReplaceScrat;
3397 ReplaceSub += SUBEXPSIZE;
3399 if ( *ma > 0 )
goto NoRep;
3400 if ( *ma <= -FUNCTION ) {
3401 *ReplaceSub++ = FUNTOFUN;
3403 *ReplaceSub++ = -*ma++;
3404 if ( *ma > -FUNCTION )
goto NoRep;
3405 *ReplaceSub++ = -*ma++;
3407 else if ( ma+4 > mb )
goto NoRep;
3409 if ( *ma == -SYMBOL ) {
3410 if ( ma[2] == -SYMBOL && ma+4 <= mb )
3411 *ReplaceSub++ = SYMTOSYM;
3412 else if ( ma[2] == -SNUMBER && ma+4 <= mb ) {
3413 *ReplaceSub++ = SYMTONUM;
3414 if ( ReplaceType == 0 ) {
3415 oldtoprhs = C->numrhs;
3420 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3421 *ReplaceSub++ = SYMTONUM;
3423 *ReplaceSub++ = ma[1];
3432 else if ( ma[2] > 0 ) {
3433 WORD *sstop, *ttstop, n;
3437 while ( ss < sstop ) {
3439 ttstop = tt - ABS(tt[-1]);
3441 while ( ss < ttstop ) {
3442 if ( *ss == INDEX )
goto NoRep;
3448 if ( ReplaceType == 0 ) {
3449 oldtoprhs = C->numrhs;
3453 ss =
AddRHS(AT.ebufnum,1);
3458 while ( --n >= 0 ) *ss++ = *tt++;
3460 C->
rhs[C->numrhs+1] = ss;
3462 *ReplaceSub++ = subtype;
3464 *ReplaceSub++ = ma[1];
3465 *ReplaceSub++ = C->numrhs;
3471 else if ( ( *ma == -VECTOR || *ma == -MINVECTOR ) && ma+4 <= mb ) {
3472 if ( ma[2] == -VECTOR ) {
3473 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOVEC;
3474 else *ReplaceSub++ = VECTOMIN;
3476 else if ( ma[2] == -MINVECTOR ) {
3477 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOMIN;
3478 else *ReplaceSub++ = VECTOVEC;
3484 else if ( ma[2] > 0 ) {
3485 WORD *sstop, *ttstop, *w, *mm, n, count;
3487 if ( *ma == -MINVECTOR ) {
3491 while ( ss < sstop ) {
3500 while ( ss < sstop ) {
3502 ttstop = tt - ABS(tt[-1]);
3505 while ( ss < ttstop ) {
3506 if ( *ss == INDEX ) {
3507 n = ss[1] - 2; ss += 2;
3508 while ( --n >= 0 ) {
3509 if ( *ss < MINSPEC ) count++;
3515 if ( count != 1 )
goto NoRep;
3519 if ( ReplaceType == 0 ) {
3520 oldtoprhs = C->numrhs;
3524 mm =
AddRHS(AT.ebufnum,1);
3525 *ReplaceSub++ = subtype;
3527 *ReplaceSub++ = ma[1];
3528 *ReplaceSub++ = C->numrhs;
3532 while ( (mm + n + 10) > C->
Top )
3534 while ( --n >= 0 ) *mm++ = *w++;
3536 C->
rhs[C->numrhs+1] = mm;
3538 mm =
AddRHS(AT.ebufnum,1);
3542 while ( (mm + n + 13) > C->
Top )
3545 while ( w < sstop ) {
3546 tt = w + *w; ttstop = tt - ABS(tt[-1]);
3548 while ( w < ttstop ) {
3549 if ( *w != INDEX ) {
3558 while ( --n >= 0 ) {
3559 if ( *w >= MINSPEC ) *mm++ = *w++;
3564 if ( n <= 2 ) mm -= 2;
3573 while ( w < tt ) *mm++ = *w++;
3574 *ss = WORDDIF(mm,ss);
3577 C->
rhs[C->numrhs+1] = mm;
3579 if ( mm > C->
Top ) {
3580 MLOCK(ErrorMessageLock);
3581 MesPrint(
"Internal error in Normalize with extra compiler buffer");
3582 MUNLOCK(ErrorMessageLock);
3590 else if ( *ma == -INDEX ) {
3591 if ( ( ma[2] == -INDEX || ma[2] == -VECTOR )
3593 *ReplaceSub++ = INDTOIND;
3594 else if ( ma[1] >= AM.OffsetIndex ) {
3595 if ( ma[2] == -SNUMBER && ma+4 <= mb
3596 && ma[3] >= 0 && ma[3] < AM.OffsetIndex )
3597 *ReplaceSub++ = INDTOIND;
3598 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3599 *ReplaceSub++ = INDTOIND;
3601 *ReplaceSub++ = ma[1];
3612 *ReplaceSub++ = ma[1];
3613 *ReplaceSub++ = ma[3];
3618 AN.ReplaceScrat[1] = ReplaceSub-AN.ReplaceScrat;
3623 while ( mb < m ) *mc++ = *mb++;
3627 if ( ReplaceType > 0 ) {
3628 C->numrhs = oldtoprhs;
3632 if ( ++ReplaceVeto >= 0 )
break;
3643 for ( i = 0; i < neps; i++ ) {
3645 if ( ( t[2] & DIRTYSYMFLAG ) != DIRTYSYMFLAG )
continue;
3646 t[2] &= ~DIRTYSYMFLAG;
3647 if ( AR.Eside == LHSIDE || AR.Eside == LHSIDEX ) {
3656 if ( *r != FUNNYWILD ) { r++;
continue; }
3657 k = r[1]; u = r + 2;
3660 if ( *u != FUNNYWILD ) ncoef = -ncoef;
3663 tt[-2] = FUNNYWILD; tt[-1] = k; m -= 2;
3667 for ( r = t + 1; r < m; r++ ) {
3668 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3669 else if ( *r == *t )
goto NormZero;
3674 for ( r = t + 2; r < tt; r += 2 ) {
3675 if ( r[1] < t[1] ) {
3676 k = r[1]; r[1] = t[1]; t[1] = k; ncoef = -ncoef; }
3677 else if ( r[1] == t[1] )
goto NormZero;
3686 for ( r = t + 1; r < m; r++ ) {
3687 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3688 else if ( *r == *t )
goto NormZero;
3697 for ( i = 0; i < (neps-1); i++ ) {
3699 for ( j = i+1; j < neps; j++ ) {
3701 if ( t[1] > r[1] ) {
3702 peps[i] = m = r; peps[j] = r = t; t = m;
3704 else if ( t[1] == r[1] ) {
3710 m = peps[j]; peps[j] = t; peps[i] = t = m;
3713 else if ( *r++ > *m++ )
break;
3714 }
while ( --k > 0 );
3719 for ( i = 0; i < neps; i++ ) {
3731 for ( i = 0; i < ndel; i += 2, r += 2 ) {
3732 if ( r[1] < r[0] ) { k = *r; *r = r[1]; r[1] = k; }
3734 for ( i = 2; i < ndel; i += 2, t += 2 ) {
3736 for ( j = i; j < ndel; j += 2 ) {
3737 if ( *r > *t ) { r += 2; }
3738 else if ( *r < *t ) {
3739 k = *r; *r++ = *t; *t++ = k;
3740 k = *r; *r++ = *t; *t-- = k;
3743 if ( *++r < t[1] ) {
3744 k = *r; *r = t[1]; t[1] = k;
3762 for ( i = 0; i < nind; i++ ) {
3764 for ( j = i+1; j < nind; j++ ) {
3766 k = *r; *r = *t; *t = k;
3784 for ( i = 2; i < nvec; i += 2 ) {
3786 for ( j = i; j < nvec; j += 2 ) {
3788 if ( *++r < t[1] ) {
3789 k = *r; *r = t[1]; t[1] = k;
3793 else if ( *r < *t ) {
3794 k = *r; *r++ = *t; *t++ = k;
3795 k = *r; *r++ = *t; *t-- = k;
3815 while ( --i >= 0 ) {
3816 if ( *t > t[1] ) { j = *t; *t = t[1]; t[1] = j; }
3822 while ( t < (m-3) ) {
3826 if ( *++r == *++t ) {
3828 if ( ( *r < MAXPOWER && t[1] < MAXPOWER )
3829 || ( *r > -MAXPOWER && t[1] > -MAXPOWER ) ) {
3832 if ( *t > MAXPOWER || *t < -MAXPOWER ) {
3833 MLOCK(ErrorMessageLock);
3834 MesPrint(
"Exponent of dotproduct out of range: %d",*t);
3835 MUNLOCK(ErrorMessageLock);
3851 else if ( *r < *++t ) {
3852 k = *r; *r++ = *t; *t = k;
3857 else if ( *r < *t ) {
3858 k = *r; *r++ = *t; *t++ = k;
3859 k = *r; *r++ = *t; *t = k;
3862 else { r += 2; t--; }
3864 else if ( *r < *t ) {
3865 k = *r; *r++ = *t; *t++ = k;
3866 k = *r; *r++ = *t; *t++ = k;
3867 k = *r; *r++ = *t; *t = k;
3876 if ( ( i = ndot ) > 0 ) {
3891 *m++ = ( i = nsym ) + 2;
3894 if ( t[1] < (2*MAXPOWER) ) {
3895 if ( t[1] & 1 ) { *m++ = 0; *m++ = 1; }
3897 if ( *++t & 2 ) ncoef = -ncoef;
3901 else if ( *t <= NumSymbols && *t > -2*MAXPOWER ) {
3902 if ( ( ( ( t[1] > symbols[*t].maxpower ) && ( symbols[*t].maxpower < MAXPOWER ) ) ||
3903 ( ( t[1] < symbols[*t].minpower ) && ( symbols[*t].minpower > -MAXPOWER ) ) ) &&
3904 ( t[1] < 2*MAXPOWER ) && ( t[1] > -2*MAXPOWER ) ) {
3905 if ( i <= 2 || t[2] != *t )
goto NormZero;
3907 if ( AN.ncmod == 1 && ( AC.modmode & ALSOPOWERS ) != 0 ) {
3908 if ( AC.cmod[0] == 1 ) t[1] = 0;
3909 else if ( t[1] >= 0 ) t[1] = 1 + (t[1]-1)%(AC.cmod[0]-1);
3911 t[1] = -1 - (-t[1]-1)%(AC.cmod[0]-1);
3912 if ( t[1] < 0 ) t[1] += (AC.cmod[0]-1);
3915 if ( ( t[1] < (2*MAXPOWER) && t[1] >= MAXPOWER )
3916 || ( t[1] > -(2*MAXPOWER) && t[1] <= -MAXPOWER ) ) {
3917 MLOCK(ErrorMessageLock);
3918 MesPrint(
"Exponent out of range: %d",t[1]);
3919 MUNLOCK(ErrorMessageLock);
3922 if ( AT.TrimPower && AR.PolyFunVar == *t && t[1] > AR.PolyFunPow ) {
3929 else { *r -= 2; t += 2; }
3932 *m++ = *t++; *m++ = *t++;
3934 }
while ( (i-=2) > 0 ); }
3935 if ( *r <= 2 ) m = r-1;
3941 stop = (WORD *)(((UBYTE *)(termout)) + AM.MaxTer);
3943 if ( ( m + i ) > stop ) {
3944 MLOCK(ErrorMessageLock);
3945 MesPrint(
"Term too complex during normalization");
3946 MUNLOCK(ErrorMessageLock);
3949 if ( ReplaceType >= 0 ) {
3956 if ( ReplaceType == 0 ) {
3957 AT.WorkPointer = termout+*termout;
3958 WildFill(BHEAD term,termout,AN.ReplaceScrat);
3959 termout = term + *term;
3962 AT.WorkPointer = r = termout + *termout;
3963 WildFill(BHEAD r,termout,AN.ReplaceScrat);
3970 r += *term; r -= ABS(r[-1]);
3973 if ( *m >= FUNCTION && m[1] > FUNHEAD &&
3974 functions[*m-FUNCTION].spec != TENSORFUNCTION )
3987 TermFree(n_llnum,
"n_llnum");
3988 TermFree(n_coef,
"NormCoef");
4007 if ( termout < term + *term && termout >= term ) AT.WorkPointer = term + *term;
4008 else AT.WorkPointer = termout;
4014 TermFree(n_llnum,
"n_llnum");
4015 TermFree(n_coef,
"NormCoef");
4019 MLOCK(ErrorMessageLock);
4020 MesPrint(
"Division by zero during normalization");
4021 MUNLOCK(ErrorMessageLock);
4025 MLOCK(ErrorMessageLock);
4026 MesPrint(
"0^0 during normalization of term");
4027 MUNLOCK(ErrorMessageLock);
4031 MLOCK(ErrorMessageLock);
4032 MesPrint(
"0/0 in polyratfun during normalization of term");
4033 MUNLOCK(ErrorMessageLock);
4038 AT.WorkPointer = termout;
4039 TermFree(n_llnum,
"n_llnum");
4040 TermFree(n_coef,
"NormCoef");
4044 TermFree(n_llnum,
"n_llnum");
4045 TermFree(n_coef,
"NormCoef");
4049 MLOCK(ErrorMessageLock);
4051 MUNLOCK(ErrorMessageLock);
4052 TermFree(n_llnum,
"n_llnum");
4053 TermFree(n_coef,
"NormCoef");
4066WORD ExtraSymbol(WORD sym, WORD pow, WORD nsym, WORD *ppsym, WORD *ncoef)
4074 if ( pow > 2*MAXPOWER || pow < -2*MAXPOWER
4075 || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
4076 MLOCK(ErrorMessageLock);
4077 MesPrint(
"Illegal wildcard power combination.");
4078 MUNLOCK(ErrorMessageLock);
4083 if ( ( sym <= NumSymbols && sym > -MAXPOWER )
4084 && ( symbols[sym].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
4085 *m %= symbols[sym].maxpower;
4086 if ( *m < 0 ) *m += symbols[sym].maxpower;
4087 if ( ( symbols[sym].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
4088 if ( ( ( symbols[sym].maxpower & 1 ) == 0 ) &&
4089 ( *m >= symbols[sym].maxpower/2 ) ) {
4090 *m -= symbols[sym].maxpower/2; *ncoef = -*ncoef;
4095 if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
4096 MLOCK(ErrorMessageLock);
4097 MesPrint(
"Power overflow during normalization");
4098 MUNLOCK(ErrorMessageLock);
4104 { *m = m[2]; m++; *m = m[2]; m++; i++; }
4109 else if ( sym < *m ) {
4117 { m--; m[2] = *m; m--; m[2] = *m; i++; }
4128WORD DoTheta(PHEAD WORD *t)
4131 WORD k, *r1, *r2, *tstop, type;
4132 WORD ia, *ta, *tb, *stopa, *stopb;
4133 if ( AC.BracketNormalize )
return(-1);
4138 if ( k <= FUNHEAD )
return(1);
4141 if ( r1 == tstop ) {
4145 if ( *t == ARGHEAD ) {
4146 if ( type == THETA )
return(1);
4150 if ( *t == -SNUMBER ) {
4151 if ( t[1] < 0 )
return(0);
4153 if ( type == THETA2 && t[1] == 0 )
return(0);
4160 if ( *t == ABS(k)+1+ARGHEAD ) {
4161 if ( k > 0 )
return(1);
4171 if ( r2 < tstop )
return(-1);
4176 if ( *t == -SNUMBER && *r1 == -SNUMBER ) {
4177 if ( t[1] > r1[1] )
return(0);
4178 else if ( t[1] < r1[1] ) {
4181 else if ( type == THETA )
return(1);
4184 else if ( t[1] == 0 && *t == -SNUMBER ) {
4186 else if ( *t < *r1 )
return(1);
4187 else if ( *t > *r1 )
return(0);
4189 else if ( r1[1] == 0 && *r1 == -SNUMBER ) {
4191 else if ( *t < *r1 )
return(1);
4192 else if ( *t > *r1 )
return(0);
4194 r2 = AT.WorkPointer;
4208 ta += ARGHEAD; tb += ARGHEAD;
4209 while ( ta < stopa ) {
4210 if ( tb >= stopb )
return(0);
4211 if ( ( ia = CompareTerms(ta,tb,(WORD)1) ) < 0 )
return(0);
4212 if ( ia > 0 )
return(1);
4216 if ( type == THETA )
return(1);
4225WORD DoDelta(WORD *t)
4227 WORD k, *r1, *r2, *tstop, isnum, isnum2, type = *t;
4228 if ( AC.BracketNormalize )
return(-1);
4230 if ( k <= FUNHEAD )
goto argzero;
4231 if ( k == FUNHEAD+ARGHEAD && t[FUNHEAD] == ARGHEAD )
goto argzero;
4238 if ( *t == -SNUMBER ) { isnum = 1; k = t[1]; }
4244 if ( k == *t-ARGHEAD-1 ) isnum = 1;
4248 if ( r1 >= tstop ) {
4249 if ( !isnum )
return(-1);
4250 if ( k == 0 )
goto argzero;
4255 if ( r2 < tstop )
return(-1);
4257 if ( *r1 == -SNUMBER ) { isnum2 = 1; }
4263 if ( k == *r1-ARGHEAD-1 ) isnum2 = 1;
4266 if ( isnum != isnum2 )
return(-1);
4268 while ( t < tstop && r1 < r2 ) {
4270 if ( !isnum )
return(-1);
4275 if ( t != tstop || r1 != r2 ) {
4276 if ( !isnum )
return(-1);
4280 if ( type == DELTA2 )
return(1);
4283 if ( type == DELTA2 )
return(0);
4292void DoRevert(WORD *fun, WORD *tmp)
4294 WORD *t, *r, *m, *to, *tt, *mm, i, j;
4299 if ( *r == -REVERSEFUNCTION ) {
4301 while ( mm < to ) *m++ = *mm++;
4304 fun[2] |= DIRTYSYMFLAG;
4306 else if ( *r <= -FUNCTION ) r++;
4308 if ( *r == -INDEX && r[1] < MINSPEC ) *r = -VECTOR;
4313 if ( ( *r > ARGHEAD )
4314 && ( r[ARGHEAD+1] == REVERSEFUNCTION )
4315 && ( *r == (r[ARGHEAD]+ARGHEAD) )
4316 && ( r[ARGHEAD] == (r[ARGHEAD+2]+4) )
4317 && ( *(r+*r-3) == 1 )
4318 && ( *(r+*r-2) == 1 )
4319 && ( *(r+*r-1) == 3 ) ) {
4331 while ( --j >= 0 ) {
4334 while ( --i >= 0 ) {
4341 else if ( *t <= -FUNCTION ) *m++ = *t++;
4342 else { *m++ = *t++; *m++ = *t++; }
4353 fun[1] = WORDDIF(t,fun);
4355 fun[2] |= DIRTYSYMFLAG;
4378#define MAXNUMBEROFNONCOMTERMS 2
4380WORD DetCommu(WORD *terms)
4382 WORD *t, *tnext, *tstop;
4384 if ( *terms == 0 )
return(0);
4385 if ( terms[*terms] == 0 )
return(0);
4389 tstop = tnext - ABS(tnext[-1]);
4391 while ( t < tstop ) {
4392 if ( *t >= FUNCTION ) {
4393 if ( functions[*t-FUNCTION].commute ) {
4395 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4399 else if ( *t == SUBEXPRESSION ) {
4400 if ( cbuf[t[4]].CanCommu[t[2]] ) {
4402 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4406 else if ( *t == EXPRESSION ) {
4408 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4411 else if ( *t == DOLLAREXPRESSION ) {
4418 if ( cbuf[AM.dbufnum].CanCommu[t[2]] ) {
4420 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4439WORD DoesCommu(WORD *term)
4443 if ( *term == 0 )
return(0);
4444 tstop = term + *term;
4445 tstop = tstop - ABS(tstop[-1]);
4447 while ( term < tstop ) {
4448 if ( ( *term >= FUNCTION ) && ( functions[*term-FUNCTION].commute ) ) {
4450 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4465WORD *PolyNormPoly (PHEAD WORD *Poly) {
4468 WORD *buffer = AT.WorkPointer;
4470 if (
NewSort(BHEAD0) ) { Terminate(-1); }
4482 if (
EndSort(BHEAD buffer,1) < 0 ) {
4487 while ( *p ) p += *p;
4489 AT.WorkPointer = p + 1;
4516WORD *EvaluateGcd(PHEAD WORD *subterm)
4519 WORD *oldworkpointer = AT.WorkPointer, *work1, *work2, *work3;
4520 WORD *t, *tt, *ttt, *t1, *t2, *t3, *t4, *tstop;
4523 WORD *lnum=n_llnum+1;
4524 WORD *num1, *num2, *num3, *den1, *den2, *den3;
4525 WORD sizenum1, sizenum2, sizenum3, sizeden1, sizeden2, sizeden3;
4526 int i, isnumeric = 0, numarg = 0 ;
4532 tt = subterm + subterm[1]; t = subterm + FUNHEAD;
4536 if ( *t == -SNUMBER ) {
4539 MLOCK(ErrorMessageLock);
4540 MesPrint(
"Trying to take the GCD involving a zero term.");
4541 MUNLOCK(ErrorMessageLock);
4545 t1 = subterm + FUNHEAD;
4546 while ( gcdnum > 1 && t1 < tt ) {
4547 if ( *t1 == -SNUMBER ) {
4549 if ( stor == 0 )
goto gcdzero;
4550 if ( GcdLong(BHEAD (UWORD *)&stor,1,(UWORD *)&gcdnum,1,
4551 (UWORD *)lnum,&nnum) )
goto FromGCD;
4556 else if ( *t1 == -SYMBOL )
goto gcdisone;
4557 else if ( *t1 < 0 )
goto gcdillegal;
4563 ct = *ttt; *ttt = 0;
4565 t1 = PolyNormPoly(BHEAD t1+ARGHEAD);
4574 while ( t3 > t2 && *t3 == 0 ) { t3--; i--; }
4575 if ( GcdLong(BHEAD (UWORD *)t2,(WORD)i,(UWORD *)&gcdnum,1,
4576 (UWORD *)lnum,&nnum) ) {
4581 if ( gcdnum == 1 ) {
4588 AT.WorkPointer = oldworkpointer;
4590 if ( gcdnum == 1 )
goto gcdisone;
4591 oldworkpointer[0] = 4;
4592 oldworkpointer[1] = gcdnum;
4593 oldworkpointer[2] = 1;
4594 oldworkpointer[3] = 3;
4595 oldworkpointer[4] = 0;
4596 AT.WorkPointer = oldworkpointer + 5;
4597 return(oldworkpointer);
4599 else if ( *t == -SYMBOL ) {
4600 t1 = subterm + FUNHEAD;
4603 if ( *t1 == -SNUMBER )
goto gcdisone;
4604 if ( *t1 == -SYMBOL ) {
4605 if ( t1[1] != i )
goto gcdisone;
4608 if ( *t1 < 0 )
goto gcdillegal;
4610 ct = *ttt; *ttt = 0;
4612 t2 = PolyNormPoly(BHEAD t1+ARGHEAD);
4614 else t2 = t1 + ARGHEAD;
4618 tstop = t2 - ABS(t2[-1]);
4619 while ( t3 < tstop ) {
4620 if ( *t3 != SYMBOL ) {
4627 if ( *t4 == i && t4[1] > 0 )
goto nextterminarg;
4637 AT.WorkPointer = oldworkpointer;
4639 oldworkpointer[0] = 8;
4640 oldworkpointer[1] = SYMBOL;
4641 oldworkpointer[2] = 4;
4642 oldworkpointer[3] = t[1];
4643 oldworkpointer[4] = 1;
4644 oldworkpointer[5] = 1;
4645 oldworkpointer[6] = 1;
4646 oldworkpointer[7] = 3;
4647 oldworkpointer[8] = 0;
4648 AT.WorkPointer = oldworkpointer+9;
4649 return(oldworkpointer);
4651 else if ( *t < 0 ) {
4653 MLOCK(ErrorMessageLock);
4654 MesPrint(
"Illegal object in gcd_ function. Object not a number or a symbol.");
4655 MUNLOCK(ErrorMessageLock);
4658 else if ( ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4659 else if ( t[1] != 0 ) {
4660 ttt = t + *t; ct = *ttt; *ttt = 0;
4661 t = PolyNormPoly(BHEAD t+ARGHEAD);
4663 if ( t[*t] == 0 && ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4664 AT.WorkPointer = oldworkpointer;
4679 AT.WorkPointer = oldworkpointer;
4681 t = subterm + FUNHEAD;
4682 for ( i = 1; i < isnumeric; i++ ) {
4686 ttt = t + *t; ct = *ttt; *ttt = 0;
4687 t = PolyNormPoly(BHEAD t+ARGHEAD);
4691 i = (ABS(t[-1])-1)/2;
4694 sizenum1 = sizeden1 = i;
4695 while ( sizenum1 > 1 && num1[sizenum1-1] == 0 ) sizenum1--;
4696 while ( sizeden1 > 1 && den1[sizeden1-1] == 0 ) sizeden1--;
4697 work1 = AT.WorkPointer+1; work2 = work1+sizenum1;
4698 for ( i = 0; i < sizenum1; i++ ) work1[i] = num1[i];
4699 for ( i = 0; i < sizeden1; i++ ) work2[i] = den1[i];
4700 num1 = work1; den1 = work2;
4701 AT.WorkPointer = work2 = work2 + sizeden1;
4702 t = subterm + FUNHEAD;
4704 ttt = t + *t; ct = *ttt; *ttt = 0;
4706 t = PolyNormPoly(BHEAD t+ARGHEAD);
4711 i = (ABS(t[-1])-1)/2;
4714 sizenum2 = sizeden2 = i;
4715 while ( sizenum2 > 1 && num2[sizenum2-1] == 0 ) sizenum2--;
4716 while ( sizeden2 > 1 && den2[sizeden2-1] == 0 ) sizeden2--;
4717 num3 = AT.WorkPointer;
4718 if ( GcdLong(BHEAD (UWORD *)num2,sizenum2,(UWORD *)num1,sizenum1,
4719 (UWORD *)num3,&sizenum3) )
goto FromGCD;
4720 sizenum1 = sizenum3;
4721 for ( i = 0; i < sizenum1; i++ ) num1[i] = num3[i];
4722 den3 = AT.WorkPointer;
4723 if ( GcdLong(BHEAD (UWORD *)den2,sizeden2,(UWORD *)den1,sizeden1,
4724 (UWORD *)den3,&sizeden3) )
goto FromGCD;
4725 sizeden1 = sizeden3;
4726 for ( i = 0; i < sizeden1; i++ ) den1[i] = den3[i];
4727 if ( sizenum1 == 1 && num1[0] == 1 && sizeden1 == 1 && den1[1] == 1 )
4732 AT.WorkPointer = work2;
4734 AT.WorkPointer = oldworkpointer;
4738 if ( sizenum1 > sizeden1 ) {
4739 while ( sizenum1 > sizeden1 ) den1[sizeden1++] = 0;
4741 else if ( sizenum1 < sizeden1 ) {
4742 while ( sizenum1 < sizeden1 ) num1[sizenum1++] = 0;
4747 if ( num1 != t ) { NCOPY(t,num1,sizenum1); }
4749 if ( den1 != t ) { NCOPY(t,den1,sizeden1); }
4754 return(oldworkpointer);
4761 t = subterm + FUNHEAD;
4762 AT.WorkPointer += AM.MaxTer/
sizeof(WORD);
4763 work2 = AT.WorkPointer;
4770 work1 = AT.WorkPointer;
4771 ttt = t + *t; ct = *ttt; *ttt = 0;
4772 t = PolyNormPoly(BHEAD t+ARGHEAD);
4773 if ( *work1 < AT.WorkPointer-work1 ) {
4782 *AT.WorkPointer++ = 0;
4788 if ( work2 != work3 ) {
4789 work1 = PolyGCD2(BHEAD work1,work2);
4791 while ( *work2 ) work2 += *work2;
4795 while ( *work2 ) work2 += *work2;
4796 size = work2 - work1 + 1;
4798 NCOPY(t,work1,size);
4800 return(oldworkpointer);
4803 oldworkpointer[0] = 4;
4804 oldworkpointer[1] = 1;
4805 oldworkpointer[2] = 1;
4806 oldworkpointer[3] = 3;
4807 oldworkpointer[4] = 0;
4808 AT.WorkPointer = oldworkpointer+5;
4809 return(oldworkpointer);
4811 MLOCK(ErrorMessageLock);
4812 MesCall(
"EvaluateGcd");
4813 MUNLOCK(ErrorMessageLock);
4830int TreatPolyRatFun(PHEAD WORD *prf)
4832 WORD *t, *tstop, *r, *rstop, *m, *mstop;
4833 WORD exp1 = MAXPOWER, exp2 = MAXPOWER;
4836 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
4837 if ( exp1 > 1 ) exp1 = 1;
4841 if ( exp1 > 0 ) exp1 = 0;
4848 while ( t < tstop ) {
4854 rstop = t - ABS(t[-1]);
4855 while ( r < rstop ) {
4856 if ( *r != SYMBOL ) { r += r[1];
continue; }
4860 while ( m < mstop ) {
4861 if ( *m == AR.PolyFunVar ) {
4862 if ( m[1] < exp1 ) exp1 = m[1];
4868 if ( exp1 > 0 ) exp1 = 0;
4873 if ( exp1 > 0 ) exp1 = 0;
4879 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
4880 if ( exp2 > 1 ) exp2 = 1;
4883 if ( exp2 > 0 ) exp2 = 0;
4889 while ( t < tstop ) {
4895 rstop = t - ABS(t[-1]);
4896 while ( r < rstop ) {
4897 if ( *r != SYMBOL ) { r += r[1];
continue; }
4901 while ( m < mstop ) {
4902 if ( *m == AR.PolyFunVar ) {
4903 if ( m[1] < exp2 ) exp2 = m[1];
4909 if ( exp2 > 0 ) exp2 = 0;
4914 if ( exp2 > 0 ) exp2 = 0;
4927 *t++ = -SNUMBER; *t++ = 1;
4928 *t++ = -SNUMBER; *t++ = 1;
4930 else if ( exp1 > 0 ) {
4932 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
4938 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
4939 *t++ = exp1; *t++ = 1; *t++ = 1; *t++ = 3;
4941 *t++ = -SNUMBER; *t++ = 1;
4944 *t++ = -SNUMBER; *t++ = 1;
4946 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
4952 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
4953 *t++ = -exp1; *t++ = 1; *t++ = 1; *t++ = 3;
4966void DropCoefficient(PHEAD WORD *term)
4969 WORD *t = term + *term;
4971 n = t[-1]; na = ABS(n);
4973 if ( n == 3 && t[0] == 1 && t[1] == 1 )
return;
4975 t[0] = 1; t[1] = 1; t[2] = 3;
4984void DropSymbols(PHEAD WORD *term)
4987 WORD *tend = term + *term, *t1, *t2, *tstop;
4988 tstop = tend - ABS(tend[-1]);
4990 while ( t1 < tstop ) {
4991 if ( *t1 == SYMBOL ) {
4994 while ( t2 < tend ) *t1++ = *t2++;
5017 WORD buffer[7*NORMSIZE], *t, *b, *bb, *tt, *m, *tstop;
5020 *b++ = SYMBOL; *b++ = 2;
5022 tstop = t - ABS(t[-1]);
5024 while ( t < tstop ) {
5025 if ( *t == SYMBOL && t < tstop ) {
5026 for ( i = 2; i < t[1]; i += 2 ) {
5029 if ( bb[0] == t[i] ) {
5031 if ( bb[1] > MAXPOWER || bb[1] < -MAXPOWER ) {
5032 MLOCK(ErrorMessageLock);
5033 MesPrint(
"Power in SymbolNormalize out of range");
5034 MUNLOCK(ErrorMessageLock);
5040 bb[0] = bb[2]; bb[1] = bb[3]; bb += 2;
5045 else if ( bb[0] > t[i] ) {
5047 while ( m > bb ) { m[1] = m[-1]; m[0] = m[-2]; m -= 2; }
5056 *b++ = t[i]; *b++ = t[i+1];
5062 MLOCK(ErrorMessageLock);
5063 MesPrint(
"Illegal term in SymbolNormalize");
5064 MUNLOCK(ErrorMessageLock);
5069 buffer[1] = b - buffer;
5073 if ( AT.LeaveNegative == 0 ) {
5074 b = buffer; bb = b + b[1]; b += 3;
5077 MLOCK(ErrorMessageLock);
5078 MesPrint(
"Negative power in SymbolNormalize");
5079 MUNLOCK(ErrorMessageLock);
5091 b = buffer; tt = term + 1;
5092 if ( i > 2 ) { NCOPY(tt,b,i) }
5095 if ( i < 0 ) i = -i;
5096 *term -= (tstop-tt);
5110int TestFunFlag(PHEAD WORD *tfun)
5112 WORD *t, *tstop, *r, *rstop, *m, *mstop;
5113 if ( functions[*tfun-FUNCTION].spec )
return(0);
5114 tstop = tfun + tfun[1];
5116 while ( t < tstop ) {
5117 if ( *t < 0 ) { NEXTARG(t);
continue; }
5119 if ( t[1] == 0 ) { t = rstop;
continue; }
5121 while ( r < rstop ) {
5122 m = r+1; mstop = r+*r; mstop -= ABS(mstop[-1]);
5123 while ( m < mstop ) {
5124 if ( *m == SUBEXPRESSION || *m == EXPRESSION || *m == DOLLAREXPRESSION )
return(1);
5125 if ( ( *m >= FUNCTION ) && ( ( m[2] & DIRTYFLAG ) == DIRTYFLAG )
5126 && ( *m != REPLACEMENT ) && TestFunFlag(BHEAD m) )
return(1);
5141WORD BracketNormalize(PHEAD WORD *term)
5143 WORD *stop = term+*term-3, *t, *tt, *tstart, *r;
5144 WORD *oldwork = AT.WorkPointer;
5147 termout = AT.WorkPointer = term+*term;
5151 tt = termout+1; t = term+1;
5152 while ( t < stop ) {
5153 if ( *t >= FUNCTION ) { i = t[1]; NCOPY(tt,t,i); }
5156 if ( tt > termout+1 && tt-termout-1 > termout[2] ) {
5157 r = termout+1; ii = tt-r;
5158 for ( i = 0; i < ii-FUNHEAD; i += FUNHEAD ) {
5159 for ( j = i+FUNHEAD; j > 0; j -= FUNHEAD ) {
5160 if ( functions[r[j-FUNHEAD]-FUNCTION].commute
5161 && functions[r[j]-FUNCTION].commute == 0 )
break;
5162 if ( r[j-FUNHEAD] > r[j] ) EXCH(r[j-FUNHEAD],r[j])
5168 tstart = tt; t = term + 1; *tt++ = DELTA; *tt++ = 2;
5169 while ( t < stop ) {
5170 if ( *t == DELTA ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5173 if ( tstart[1] > 2 ) {
5174 for ( r = tstart+2; r < tstart+tstart[1]; r += 2 ) {
5175 if ( r[0] > r[1] ) EXCH(r[0],r[1])
5178 if ( tstart[1] > 4 ) {
5179 r = tstart+2; ii = tstart[1]-2;
5180 for ( i = 0; i < ii-2; i += 2 ) {
5181 for ( j = i+2; j > 0; j -= 2 ) {
5182 if ( r[j-2] > r[j] ) {
5186 else if ( r[j-2] < r[j] )
break;
5188 if ( r[j-1] > r[j+1] ) EXCH(r[j-1],r[j+1])
5193 tt = tstart+tstart[1];
5195 else if ( tstart[1] == 2 ) { tt = tstart; }
5198 tstart = tt; t = term + 1; *tt++ = INDEX; *tt++ = 2;
5199 while ( t < stop ) {
5200 if ( *t == INDEX ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5203 if ( tstart[1] >= 4 ) {
5204 r = tstart+2; ii = tstart[1]-2;
5205 for ( i = 0; i < ii-1; i += 1 ) {
5206 for ( j = i+1; j > 0; j -= 1 ) {
5207 if ( r[j-1] > r[j] ) EXCH(r[j-1],r[j])
5211 tt = tstart+tstart[1];
5213 else if ( tstart[1] == 2 ) { tt = tstart; }
5216 tstart = tt; t = term + 1; *tt++ = DOTPRODUCT; *tt++ = 2;
5217 while ( t < stop ) {
5218 if ( *t == DOTPRODUCT ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5221 if ( tstart[1] > 5 ) {
5222 r = tstart+2; ii = tstart[1]-2;
5223 for ( i = 0; i < ii; i += 3 ) {
5224 if ( r[i] > r[i+1] ) EXCH(r[i],r[i+1])
5226 for ( i = 0; i < ii-3; i += 3 ) {
5227 for ( j = i+3; j > 0; j -= 3 ) {
5228 if ( r[j-3] < r[j] )
break;
5229 if ( r[j-3] > r[j] ) {
5234 if ( r[j-2] > r[j+1] ) EXCH(r[j-2],r[j+1])
5239 tt = tstart+tstart[1];
5241 else if ( tstart[1] == 2 ) { tt = tstart; }
5243 if ( tstart[2] > tstart[3] ) EXCH(tstart[2],tstart[3])
5247 tstart = tt; t = term + 1; *tt++ = SYMBOL; *tt++ = 2;
5248 while ( t < stop ) {
5249 if ( *t == SYMBOL ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5252 if ( tstart[1] > 4 ) {
5253 r = tstart+2; ii = tstart[1]-2;
5254 for ( i = 0; i < ii-2; i += 2 ) {
5255 for ( j = i+2; j > 0; j -= 2 ) {
5256 if ( r[j-2] > r[j] ) EXCH(r[j-2],r[j])
5260 tt = tstart+tstart[1];
5262 else if ( tstart[1] == 2 ) { tt = tstart; }
5265 tstart = tt; t = term + 1; *tt++ = SETSET; *tt++ = 2;
5266 while ( t < stop ) {
5267 if ( *t == SETSET ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5270 if ( tstart[1] > 4 ) {
5271 r = tstart+2; ii = tstart[1]-2;
5272 for ( i = 0; i < ii-2; i += 2 ) {
5273 for ( j = i+2; j > 0; j -= 2 ) {
5274 if ( r[j-2] > r[j] ) {
5281 tt = tstart+tstart[1];
5283 else if ( tstart[1] == 2 ) { tt = tstart; }
5285 *tt++ = 1; *tt++ = 1; *tt++ = 3;
5286 t = term; i = *termout = tt - termout; tt = termout;
5288 AT.WorkPointer = oldwork;
WORD CompareSymbols(WORD *, WORD *, WORD)
WORD CompCoef(WORD *, WORD *)
LONG EndSort(PHEAD WORD *, int)
WORD StoreTerm(PHEAD WORD *)
WORD NextPrime(PHEAD WORD)
WORD Compare1(WORD *, WORD *, WORD)
int SymbolNormalize(WORD *term)