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, 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]);
1135 multermnum:
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++; }
1339 signogood: 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]) ) {
1367 absnosymbols: 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;
1393 absnogood: 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 ) {
1737 gcdcalc:
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);
1931 nospec: 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;
1963 ncoef = REDLENG(ncoef);
1964 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(t+3),t[2]) )
goto FromNorm;
1965 ncoef = INCLENG(ncoef);
1970 if ( t[3] & 1 ) ncoef = -ncoef;
1972 else if ( t[2] == 0 ) {
1973 if ( t[3] < 0 )
goto NormInf;
1978 if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) )
goto FromNorm;
1979 ncoef = REDLENG(ncoef);
1981 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1984 else if ( t[3] > 0 ) {
1985 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1988 ncoef = INCLENG(ncoef);
1995 if ( t[1] == FUNHEAD ) {
1996 MLOCK(ErrorMessageLock);
1997 MesPrint(
"Gamma matrix without spin line encountered.");
1998 MUNLOCK(ErrorMessageLock);
2006 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG ) {
2008 t[2] |= DIRTYSYMFLAG;
2011 ScanCont:
while ( t < r ) {
2012 if ( *t >= AM.OffsetIndex &&
2013 ( *t >= AM.DumInd || ( *t < AM.WilInd &&
2014 indices[*t-AM.OffsetIndex].dimension ) ) )
2024 if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) )
2026 if ( *rr == -SNUMBER && rr[1] == 1 )
break;
2027 if ( *rr <= -FUNCTION ) k = *rr;
2029 if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) ) {
2030 if ( k == 0 )
goto NormZZ;
2033 if ( *rr == -SNUMBER && rr[1] > 0 && rr[1] < MAXPOWER && k < 0 ) {
2035 if ( functions[k-FUNCTION].commute ) {
2036 for ( i = 0; i < rr[1]; i++ ) pnco[nnco++] = rr-1;
2039 for ( i = 0; i < rr[1]; i++ ) pcom[ncom++] = rr-1;
2043 if ( k == 0 )
goto NormZero;
2044 if ( t[FUNHEAD] == -SYMBOL && *rr == -SNUMBER && t[1] == FUNHEAD+4 ) {
2045 if ( rr[1] < MAXPOWER ) {
2046 t[FUNHEAD+2] = t[FUNHEAD+1]; t += FUNHEAD+2;
2056 t[2] &= ~DIRTYSYMFLAG;
2062 t[2] &= ~DIRTYSYMFLAG;
2069 if ( *t == 0 || *t == AM.vectorzero )
goto NormZero;
2070 if ( *t > 0 && *t < AM.OffsetIndex ) {
2073 ncoef = REDLENG(ncoef);
2074 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2076 ncoef = INCLENG(ncoef);
2078 else if ( *t == NOINDEX ) t++;
2079 else pind[nind++] = *t++;
2082 case SUBEXPRESSION :
2083 if ( t[3] == 0 )
break;
2095 if ( t[2] == 0 )
goto defaultcase;
2096 if ( t[FUNHEAD] != -SNUMBER || t[FUNHEAD+1] < 0 )
goto defaultcase;
2097 if ( t[FUNHEAD+2] == -SNUMBER ) {
2098 if ( t[FUNHEAD+1] == 0 && t[FUNHEAD+3] == 0 )
goto NormZZ;
2099 if ( t[FUNHEAD+1] == 0 )
break;
2100 if ( t[FUNHEAD+3] < 0 ) {
2101 AT.WorkPointer[0] = -t[FUNHEAD+3];
2105 AT.WorkPointer[0] = t[FUNHEAD+3];
2108 AT.WorkPointer[1] = 1;
2110 else if ( t[FUNHEAD+2] == t[1]-FUNHEAD-2
2111 && t[FUNHEAD+2] == t[FUNHEAD+2+ARGHEAD]+ARGHEAD
2112 && ABS(t[t[1]-1]) == t[FUNHEAD+2+ARGHEAD] - 1 ) {
2114 if ( t[FUNHEAD+1] == 0 )
break;
2115 i = t[t[1]-1]; r1 = t + FUNHEAD+ARGHEAD+3;
2118 r2 = AT.WorkPointer;
2119 while ( --i >= 0 ) *r2++ = *r1++;
2121 else goto defaultcase;
2122 if ( TakeRatRoot((UWORD *)AT.WorkPointer,&nc,t[FUNHEAD+1]) ) {
2126 ncoef = REDLENG(ncoef);
2127 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,nc) )
2129 if ( nc < 0 ) nc = -nc;
2130 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(AT.WorkPointer+nc),nc) )
2132 ncoef = INCLENG(ncoef);
2135 case RANDOMFUNCTION :
2137 WORD nnc, nc, nca, nr;
2148 if ( t[1] == FUNHEAD )
goto defaultcase;
2149 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER &&
2150 t[FUNHEAD+1] > 0 ) {
2151 if ( t[FUNHEAD+1] == 1 )
break;
2153 ((UWORD *)AT.WorkPointer)[0] = wranf(BHEAD0);
2154 ((UWORD *)AT.WorkPointer)[1] = wranf(BHEAD0);
2156 if ( ((UWORD *)AT.WorkPointer)[1] == 0 ) {
2158 if ( ((UWORD *)AT.WorkPointer)[0] == 0 ) {
2162 xx = (UWORD)(t[FUNHEAD+1]);
2164 DivLong((UWORD *)AT.WorkPointer,nr
2166 ,((UWORD *)AT.WorkPointer)+4,&nnc
2167 ,((UWORD *)AT.WorkPointer)+2,&nc);
2168 ((UWORD *)AT.WorkPointer)[4] = 0;
2169 ((UWORD *)AT.WorkPointer)[5] = 0;
2170 ((UWORD *)AT.WorkPointer)[6] = 1;
2171 DivLong((UWORD *)AT.WorkPointer+4,3
2173 ,((UWORD *)AT.WorkPointer)+9,&nnc
2174 ,((UWORD *)AT.WorkPointer)+7,&nca);
2175 AddLong((UWORD *)AT.WorkPointer+4,3
2176 ,((UWORD *)AT.WorkPointer)+7,-nca
2177 ,((UWORD *)AT.WorkPointer)+9,&nnc);
2178 if ( BigLong((UWORD *)AT.WorkPointer,nr
2179 ,((UWORD *)AT.WorkPointer)+9,nnc) >= 0 )
goto redoshort;
2183 AT.WorkPointer[2] = (WORD)xx;
2186 ncoef = REDLENG(ncoef);
2187 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+2,nc) )
2189 ncoef = INCLENG(ncoef);
2191 else if ( t[FUNHEAD] > 0 && t[1] == t[FUNHEAD]+FUNHEAD
2192 && ABS(t[t[1]-1]) == t[FUNHEAD]-1-ARGHEAD && t[t[1]-1] > 0 ) {
2193 WORD nna, nnb, nni, nnb2, nnb2a;
2198 nnt = (UWORD *)(t+t[1]-1-nnb);
2199 if ( *nnt != 1 )
goto defaultcase;
2200 for ( nni = 1; nni < nnb; nni++ ) {
2201 if ( nnt[nni] != 0 )
goto defaultcase;
2203 nnt = (UWORD *)(t + FUNHEAD + ARGHEAD + 1);
2205 for ( nni = 0; nni < nnb2; nni++ ) {
2206 ((UWORD *)AT.WorkPointer)[nni] = wranf(BHEAD0);
2209 while ( nnb2a > 0 && ((UWORD *)AT.WorkPointer)[nnb2a-1] == 0 ) nnb2a--;
2211 DivLong((UWORD *)AT.WorkPointer,nnb2a
2213 ,((UWORD *)AT.WorkPointer)+2*nnb2,&nnc
2214 ,((UWORD *)AT.WorkPointer)+nnb2,&nc);
2215 for ( nni = 0; nni < nnb2; nni++ ) {
2216 ((UWORD *)AT.WorkPointer)[nni+2*nnb2] = 0;
2218 ((UWORD *)AT.WorkPointer)[3*nnb2] = 1;
2219 DivLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2221 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc
2222 ,((UWORD *)AT.WorkPointer)+3*nnb2+1,&nca);
2223 AddLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2224 ,((UWORD *)AT.WorkPointer)+3*nnb2+1,-nca
2225 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc);
2226 if ( BigLong((UWORD *)AT.WorkPointer,nnb2a
2227 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,nnc) >= 0 )
goto redoshort;
2231 for ( nni = 0; nni < nnb; nni++ ) {
2232 ((UWORD *)AT.WorkPointer)[nnb2+nni] = nnt[nni];
2236 ncoef = REDLENG(ncoef);
2237 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+nnb2,nc) )
2239 ncoef = INCLENG(ncoef);
2241 else goto defaultcase;
2245 if ( *t == RANPERM && t[1] > FUNHEAD && t[FUNHEAD] <= -FUNCTION ) {
2247 WORD *mm, *ww, *ow = AT.WorkPointer;
2248 WORD *Array, *targ, *argstop, narg = 0, itot;
2252 while ( targ < argstop ) {
2253 narg++; NEXTARG(targ);
2255 WantAddPointers(narg);
2256 pwork = AT.pWorkSpace + AT.pWorkPointer;
2257 targ = t+FUNHEAD+1; narg = 0;
2258 while ( targ < argstop ) {
2259 pwork[narg++] = targ;
2266 ow = AT.WorkPointer;
2267 Array = AT.WorkPointer;
2268 AT.WorkPointer += narg;
2269 for ( i = 0; i < narg; i++ ) Array[i] = i;
2270 for ( i = 2; i <= narg; i++ ) {
2271 itot = (WORD)(iranf(BHEAD i));
2272 for ( j = 0; j < itot; j++ ) CYCLE1(WORD,Array,i)
2274 mm = AT.WorkPointer;
2275 *mm++ = -t[FUNHEAD];
2277 for ( ie = 2; ie < FUNHEAD; ie++ ) *mm++ = t[ie];
2278 for ( i = 0; i < narg; i++ ) {
2279 ww = pwork[Array[i]];
2282 mm = AT.WorkPointer; t++; ww = t;
2283 i = mm[1]; NCOPY(ww,mm,i)
2284 AT.WorkPointer = ow;
2290 if ( ( t[2] & DIRTYFLAG ) != 0 && t[FUNHEAD] <= -FUNCTION
2291 && t[FUNHEAD+1] == -SNUMBER && t[FUNHEAD+2] > 0 ) {
2292 WORD *rr = t+t[1], *mm = t+FUNHEAD+3, *tt, *tt1, *tt2, num = 0;
2296 while ( mm < rr ) { num++; NEXTARG(mm); }
2297 if ( num < t[FUNHEAD+2] ) { pnco[nnco++] = t;
break; }
2298 *t = -t[FUNHEAD]; mm = t+FUNHEAD+3;
2300 while ( --i > 0 ) { NEXTARG(mm); }
2301 tt = TermMalloc(
"Select_");
2304 for ( i = 0; i < *mm; i++ ) *tt1++ = mm[i];
2306 else if ( *mm <= -FUNCTION ) { *tt1++ = *mm; }
2307 else { *tt1++ = mm[0]; *tt1++ = mm[1]; }
2309 while ( tt2 < mm ) *tt1++ = *tt2++;
2310 i = tt1-tt; tt1 = tt; tt2 = t+FUNHEAD;
2312 TermFree(tt,
"Select_");
2314 while ( mm < rr ) *tt2++ = *mm++;
2317 while ( mm < rr ) *tt2++ = *mm++;
2321 else pnco[nnco++] = t;
2332 if ( t[1] <= FUNHEAD )
break;
2337 if ( *to == ARGHEAD )
goto NormZero;
2341 if ( to[ARGHEAD] != j+1 )
goto NoInteg;
2342 if ( rr >= r ) k = -1;
2343 else if ( *rr == ARGHEAD ) { k = 0; rr += ARGHEAD; }
2344 else if ( *rr == -SNUMBER ) { k = rr[1]; rr += 2; }
2346 if ( rr != r )
goto NoInteg;
2347 if ( k > 1 || k < -1 )
goto NoInteg;
2350 i = ( i < 0 ) ? -j: j;
2351 UnPack((UWORD *)to,i,&den,&num);
2359 if ( AN.NoScrat2 == 0 ) {
2360 AN.NoScrat2 = (UWORD *)Malloc1((AM.MaxTal+2)*
sizeof(UWORD),
"Normalize");
2362 if ( DivLong((UWORD *)to,num,(UWORD *)(to+j),den
2363 ,(UWORD *)AT.WorkPointer,&num,AN.NoScrat2,&den) )
goto FromNorm;
2364 if ( k < 0 && den < 0 ) {
2367 if ( AddLong((UWORD *)AT.WorkPointer,num
2368 ,AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2371 else if ( k > 0 && den > 0 ) {
2374 if ( AddLong((UWORD *)AT.WorkPointer,num,
2375 AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2380 else if ( *to == -SNUMBER ) {
2381 if ( to[1] < 0 ) { *AT.WorkPointer = -to[1]; num = -1; }
2382 else if ( to[1] == 0 )
goto NormZero;
2383 else { *AT.WorkPointer = to[1]; num = 1; }
2386 if ( num == 0 )
goto NormZero;
2387 ncoef = REDLENG(ncoef);
2388 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,num) )
2390 ncoef = INCLENG(ncoef);
2399 if ( *t < FUNCTION ) {
2400 MLOCK(ErrorMessageLock);
2401 MesPrint(
"Illegal code in Norm");
2405 AO.OutFill = AO.OutputLine = OutBuf;
2410 while ( --i >= 0 ) {
2411 TalToLine((UWORD)(*t++));
2412 TokenToLine((UBYTE *)
" ");
2418 MUNLOCK(ErrorMessageLock);
2421 if ( *t == REPLACEMENT ) {
2422 if ( AR.Eside != LHSIDE ) ReplaceVeto--;
2430 if ( *t == DUMMYFUN || *t == DUMMYTEN ) {}
2432 if ( *t < (FUNCTION + WILDOFFSET) ) {
2433 if ( ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2434 || ( functions[*t-FUNCTION].minnumargs > 0 ) )
2435 && ( ( t[2] & DIRTYFLAG ) != 0 ) ) {
2439 WORD *ta = t + FUNHEAD, *tb = t + t[1];
2441 while ( ta < tb ) { numarg++; NEXTARG(ta) }
2442 if ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2443 && ( numarg >= functions[*t-FUNCTION].maxnumargs ) )
2445 if ( ( functions[*t-FUNCTION].minnumargs > 0 )
2446 && ( numarg < functions[*t-FUNCTION].minnumargs ) )
2450 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION].tabl == 0 ) ) {
2452 t[2] |= DIRTYSYMFLAG;
2454 if ( functions[*t-FUNCTION].commute ) { pnco[nnco++] = t; }
2455 else { pcom[ncom++] = t; }
2458 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION-WILDOFFSET].tabl == 0 ) ) {
2460 t[2] |= DIRTYSYMFLAG;
2462 if ( functions[*t-FUNCTION-WILDOFFSET].commute ) {
2465 else { pcom[ncom++] = t; }
2471 if ( ( *t < (FUNCTION + WILDOFFSET)
2472 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) || (
2473 *t >= (FUNCTION + WILDOFFSET)
2474 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION ) ) {
2475 if ( *t >= GAMMA && *t <= GAMMASEVEN ) t++;
2478 if ( *t == AM.vectorzero )
goto NormZero;
2479 if ( *t >= AM.OffsetIndex && ( *t >= AM.DumInd
2480 || ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2483 else if ( *t == FUNNYWILD ) { t++; }
2498 else if ( *t <= -FUNCTION ) t++;
2499 else if ( *t == -INDEX ) {
2500 if ( t[1] >= AM.OffsetIndex &&
2501 ( t[1] >= AM.DumInd || ( t[1] < AM.WilInd
2502 && indices[t[1]-AM.OffsetIndex].dimension ) ) )
2506 else if ( *t == -SYMBOL ) {
2507 if ( t[1] >= MAXPOWER && t[1] < 2*MAXPOWER ) {
2511 else if ( t[1] < -MAXPOWER && t[1] > -2*MAXPOWER ) {
2527 r = t = ANsr; m = ANsm;
2528 ANsc = ANsm = ANsr = 0;
2541 for ( k = 0, i = 0; i < nden; i++ ) {
2543 if ( ( t[2] & DIRTYFLAG ) == 0 )
continue;
2544 r = t + t[1]; m = t + FUNHEAD;
2546 for ( j = i+1; j < nden; j++ ) pden[j-1] = pden[j];
2548 for ( j = 0; j < nnco; j++ )
if ( pnco[j] == t )
break;
2549 for ( j++; j < nnco; j++ ) pnco[j-1] = pnco[j];
2555 if ( m >= r )
continue;
2560 k = 1; to = termout; from = term;
2562 while ( from < t ) *to++ = *from++;
2566 *to++ = DENOMINATOR;
2567 for ( j = 1; j < FUNHEAD; j++ ) *to++ = 0;
2568 if ( *m < -FUNCTION ) *to++ = *m++;
2569 else if ( *m < 0 ) { *to++ = *m++; *to++ = *m++; }
2571 j = *m;
while ( --j >= 0 ) *to++ = *m++;
2573 stop[1] = WORDDIF(to,stop);
2576 if ( i == nden - 1 ) {
2577 stop = term + *term;
2578 while ( from < stop ) *to++ = *from++;
2579 i = *termout = WORDDIF(to,termout);
2580 to = term; from = termout;
2581 while ( --i >= 0 ) *to++ = *from++;
2586 for ( i = 0; i < nden; i++ ) {
2588 if ( ( t[2] & DIRTYFLAG ) == 0 )
continue;
2590 if ( t[FUNHEAD] == -SYMBOL ) {
2593 change = ExtraSymbol(*t,-1,nsym,ppsym,&ncoef);
2595 ppsym += change * 2;
2598 else if ( t[FUNHEAD] == -SNUMBER ) {
2600 if ( *t == 0 )
goto NormInf;
2601 if ( *t < 0 ) { *AT.WorkPointer = -*t; j = -1; }
2602 else { *AT.WorkPointer = *t; j = 1; }
2603 ncoef = REDLENG(ncoef);
2604 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,j) )
2606 ncoef = INCLENG(ncoef);
2609 else if ( t[FUNHEAD] == ARGHEAD )
goto NormInf;
2610 else if ( t[FUNHEAD] > 0 && t[FUNHEAD+ARGHEAD] ==
2611 t[FUNHEAD]-ARGHEAD ) {
2614 t += FUNHEAD + ARGHEAD + 1;
2616 m = r - ABS(*r) + 1;
2617 if ( j != 3 || ( ( *m != 1 ) || ( m[1] != 1 ) ) ) {
2618 ncoef = REDLENG(ncoef);
2619 if ( DivRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,REDLENG(j),(UWORD *)n_coef,&ncoef) )
goto FromNorm;
2620 ncoef = INCLENG(ncoef);
2622 t[-FUNHEAD-ARGHEAD] -= j;
2630 if ( *t == SYMBOL || *t == DOTPRODUCT ) {
2633 pden[i][FUNHEAD] -= k;
2634 pden[i][FUNHEAD+ARGHEAD] -= k;
2639 if ( *t == SYMBOL ) {
2643 change = ExtraSymbol(*t,-t[1],nsym,ppsym,&ncoef);
2645 ppsym += change * 2;
2658 while ( to < stop ) *to++ = *from++;
2663 if ( pden[i][1] == 4+FUNHEAD+ARGHEAD ) {
2665 for ( j = 0; j < nnco; j++ ) {
2666 if ( pden[i] == pnco[j] ) {
2668 while ( j < nnco ) {
2669 pnco[j] = pnco[j+1];
2675 pden[i--] = pden[--nden];
2686 for ( i = 0; i < ndel; i += 2 ) {
2687 if ( t[0] == t[1] ) {
2688 if ( t[0] == EMPTYINDEX ) {}
2689 else if ( *t < AM.OffsetIndex ) {
2690 k = AC.FixIndices[*t];
2691 if ( k < 0 ) { j = -1; k = -k; }
2692 else if ( k > 0 ) j = 1;
2696 else if ( *t >= AM.DumInd ) {
2698 if ( k )
goto docontract;
2700 else if ( *t >= AM.WilInd ) {
2701 k = indices[*t-AM.OffsetIndex-WILDOFFSET].dimension;
2702 if ( k )
goto docontract;
2704 else if ( ( k = indices[*t-AM.OffsetIndex].dimension ) != 0 ) {
2708 WithFix: shortnum = k;
2709 ncoef = REDLENG(ncoef);
2710 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),j) )
2712 ncoef = INCLENG(ncoef);
2716 change = ExtraSymbol((WORD)(-k),(WORD)1,nsym,ppsym,&ncoef);
2718 ppsym += change * 2;
2720 t[1] = pdel[ndel-1];
2721 t[0] = pdel[ndel-2];
2728 if ( *t < AM.OffsetIndex && t[1] < AM.OffsetIndex )
goto NormZero;
2729 j = *t - AM.OffsetIndex;
2730 if ( j >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2731 || ( *t < AM.WilInd && indices[j].dimension ) ) ) {
2732 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2735 *m++ = pdel[ndel-2];
2739 else if ( *t == m[1] ) {
2741 *m++ = pdel[ndel-2];
2747 j = t[1]-AM.OffsetIndex;
2748 if ( j >= 0 && ( ( t[1] >= AM.DumInd && AC.lDefDim )
2749 || ( t[1] < AM.WilInd && indices[j].dimension ) ) ) {
2750 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2753 *m++ = pdel[ndel-2];
2757 else if ( t[1] == m[1] ) {
2759 *m++ = pdel[ndel-2];
2771 for ( i = 0; i < ndel; i++ ) {
2772 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2773 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2775 for ( j = 1; j < nvec; j += 2 ) {
2779 *t-- = pdel[--ndel];
2784 t[1] = pdel[--ndel];
2787 *t-- = pdel[--ndel];
2796 if ( ndel > 0 && ncon ) {
2798 for ( i = 0; i < ndel; i++ ) {
2799 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2800 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2801 for ( j = 0; j < ncon; j++ ) {
2802 if ( *pcon[j] == *t ) {
2805 *t-- = pdel[--ndel];
2810 t[1] = pdel[--ndel];
2813 *t-- = pdel[--ndel];
2816 for ( j = 0; j < nnco; j++ ) {
2818 if ( r > m && r < m+m[1] ) {
2819 m[2] |= DIRTYSYMFLAG;
2823 for ( j = 0; j < ncom; j++ ) {
2825 if ( r > m && r < m+m[1] ) {
2826 m[2] |= DIRTYSYMFLAG;
2830 for ( j = 0; j < neps; j++ ) {
2832 if ( r > m && r < m+m[1] ) {
2833 m[2] |= DIRTYSYMFLAG;
2848 for ( i = 3; i < nvec; i += 2 ) {
2849 k = *t - AM.OffsetIndex;
2850 if ( k >= 0 && ( ( *t > AM.DumInd && AC.lDefDim )
2851 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2853 for ( j = i; j < nvec; j += 2 ) {
2859 *r-- = pvec[--nvec];
2861 *t-- = pvec[--nvec];
2862 *t-- = pvec[--nvec];
2871 if ( nvec > 0 && ncon ) {
2873 for ( i = 1; i < nvec; i += 2 ) {
2874 k = *t - AM.OffsetIndex;
2875 if ( k >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2876 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2877 for ( j = 0; j < ncon; j++ ) {
2878 if ( *pcon[j] == *t ) {
2880 *t-- = pvec[--nvec];
2881 *t-- = pvec[--nvec];
2883 pcon[j] = pcon[--ncon];
2885 for ( j = 0; j < nnco; j++ ) {
2887 if ( r > m && r < m+m[1] ) {
2888 m[2] |= DIRTYSYMFLAG;
2892 for ( j = 0; j < ncom; j++ ) {
2894 if ( r > m && r < m+m[1] ) {
2895 m[2] |= DIRTYSYMFLAG;
2899 for ( j = 0; j < neps; j++ ) {
2901 if ( r > m && r < m+m[1] ) {
2902 m[2] |= DIRTYSYMFLAG;
2920 for ( i = 0; i < nnco; i++ ) {
2922 if ( ( *t >= (FUNCTION+WILDOFFSET)
2923 && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
2924 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2925 && functions[*t-FUNCTION].spec == 0 ) ) {
2931 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
2934 pnco[i][2] |= DIRTYSYMFLAG;
2946 for ( i = 0; i < nnco; i++ ) {
2948 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) && *t != DOLLAREXPRESSION ) {
2950 if ( ( *t >= (FUNCTION+WILDOFFSET)
2951 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
2952 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2953 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
2954 if ( *t >= (FUNCTION+WILDOFFSET) ) {
2956 j = FullSymmetrize(BHEAD t,l);
2959 else j = FullSymmetrize(BHEAD t,l);
2960 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
2961 if ( ( j & 2 ) != 0 )
goto NormZero;
2962 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
2965 else t[2] &= ~DIRTYSYMFLAG;
2973 for ( i = 0; i < k; i++ ) {
2975 while ( Commute(pnco[j],pnco[j+1]) ) {
2976 t = pnco[j]; pnco[j] = pnco[j+1]; pnco[j+1] = t;
2978 while ( l >= 0 && Commute(pnco[l],pnco[l+1]) ) {
2979 t = pnco[l]; pnco[l] = pnco[l+1]; pnco[l+1] = t;
2982 if ( ++j >= k )
break;
2989 for ( i = 0; i < nnco; i++ ) {
2991 if ( *t == IDFUNCTION ) AN.idfunctionflag = 1;
2992 if ( *t >= GAMMA && *t <= GAMMASEVEN ) {
2998 *m++ = stype = t[FUNHEAD];
3003 if ( *t == GAMMAFIVE ) {
3004 gtype = GAMMA5; t += FUNHEAD;
goto onegammamatrix; }
3005 else if ( *t == GAMMASIX ) {
3006 gtype = GAMMA6; t += FUNHEAD;
goto onegammamatrix; }
3007 else if ( *t == GAMMASEVEN ) {
3008 gtype = GAMMA7; t += FUNHEAD;
goto onegammamatrix; }
3013 if ( gtype == GAMMA5 ) {
3014 if ( j == GAMMA1 ) j = GAMMA5;
3015 else if ( j == GAMMA5 ) j = GAMMA1;
3016 else if ( j == GAMMA7 ) ncoef = -ncoef;
3017 if ( nnum & 1 ) ncoef = -ncoef;
3019 else if ( gtype == GAMMA6 || gtype == GAMMA7 ) {
3021 if ( gtype == GAMMA6 ) gtype = GAMMA7;
3022 else gtype = GAMMA6;
3024 if ( j == GAMMA1 ) j = gtype;
3025 else if ( j == GAMMA5 ) {
3027 if ( j == GAMMA7 ) ncoef = -ncoef;
3029 else if ( j != gtype )
goto NormZero;
3032 ncoef = REDLENG(ncoef);
3033 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),1) )
goto FromNorm;
3034 ncoef = INCLENG(ncoef);
3038 *m++ = gtype; nnum++;
3043 }
while ( ( ++i < nnco ) && ( *(t = pnco[i]) >= GAMMA
3044 && *t <= GAMMASEVEN ) && ( t[FUNHEAD] == stype ) );
3047 k = WORDDIF(m,to) - FUNHEAD-1;
3050 while ( --k >= 0 ) *from-- = *--r;
3053 to[1] = WORDDIF(m,to);
3055 else if ( *t < 0 ) {
3056 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3060 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3061 && *t != REPLACEMENT && *t != DOLLAREXPRESSION
3062 && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3074 for ( i = 0; i < ncom; i++ ) {
3076 if ( ( *t >= (FUNCTION+WILDOFFSET)
3077 && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
3078 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3079 && functions[*t-FUNCTION].spec == 0 ) ) {
3085 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3088 pcom[i][2] |= DIRTYSYMFLAG;
3098 for ( i = 0; i < ncom; i++ ) {
3100 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) ) {
3102 if ( ( *t >= (FUNCTION+WILDOFFSET)
3103 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3104 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3105 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3106 if ( *t >= (FUNCTION+WILDOFFSET) ) {
3108 j = FullSymmetrize(BHEAD t,l);
3111 else j = FullSymmetrize(BHEAD t,l);
3112 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3113 if ( ( j & 2 ) != 0 )
goto NormZero;
3114 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3117 else t[2] &= ~DIRTYSYMFLAG;
3126 for ( i = 1; i < ncom; i++ ) {
3127 for ( j = i; j > 0; j-- ) {
3133 if ( *r < 0 ) {
if ( *t >= *r )
goto NextI; }
3134 else {
if ( -*t <= *r )
goto NextI; }
3137 else if ( *r < 0 ) {
3138 if ( *t < -*r )
goto NextI;
3141 else if ( *t != *r ) {
3142 if ( *t < *r )
goto NextI;
3143 jexch: t = pcom[j]; pcom[j] = pcom[jj]; pcom[jj] = t;
3146 if ( AC.properorderflag ) {
3147 if ( ( *t >= (FUNCTION+WILDOFFSET)
3148 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
3149 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3150 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) ) {}
3152 WORD *s1, *s2, *ss1, *ss2;
3153 s1 = t+FUNHEAD; s2 = r+FUNHEAD;
3154 ss1 = t + t[1]; ss2 = r + r[1];
3155 while ( s1 < ss1 && s2 < ss2 ) {
3157 if ( k > 0 )
goto jexch;
3158 if ( k < 0 )
goto NextI;
3162 if ( s1 < ss1 ) goto jexch;
3166 kk = r[1] - FUNHEAD;
3169 while ( k > 0 && kk > 0 ) {
3170 if ( *t < *r )
goto NextI;
3171 else if ( *t++ > *r++ )
goto jexch;
3174 if ( k > 0 )
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 for ( i = 0; i < ncom; i++ ) {
3196 if ( *t == THETA || *t == THETA2 ) {
3197 if ( ( k = DoTheta(BHEAD t) ) == 0 )
goto NormZero;
3203 else if ( *t == DELTA2 || *t == DELTAP ) {
3204 if ( ( k = DoDelta(t) ) == 0 )
goto NormZero;
3210 else if ( *t == AR.PolyFunInv && AR.PolyFunType == 2 ) {
3215 WORD *mm, *tt = t, numt = 0;
3217 while ( tt < t+t[1] ) { numt++; NEXTARG(tt) }
3219 tt = t; mm = m; k = t[1];
3225 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3226 else { *tt++ = *mm++; *tt++ = *mm++; }
3229 k = *mm; NCOPY(tt,mm,k)
3233 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3234 else { *tt++ = *mm++; *tt++ = *mm++; }
3237 k = *mm; NCOPY(tt,mm,k)
3240 t[2] |= MUSTCLEANPRF;
3244 else if ( *t == AR.PolyFun ) {
3245 if ( AR.PolyFunType == 1 ) {
3246 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3247 t[1] == FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER )
goto NormZero;
3248 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) {
3249 if ( AN.PolyNormFlag == 0 ) {
3250 AN.PolyNormFlag = 1;
3257 else if ( AR.PolyFunType == 2 ) {
3268 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3269 t[1] > FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) {
3270 u = t + FUNHEAD + 2;
3272 if ( *u <= -FUNCTION ) {}
3273 else if ( t[1] == FUNHEAD+4 && t[FUNHEAD+2] == -SNUMBER
3274 && t[FUNHEAD+3] == 0 )
goto NormPRF;
3275 else if ( t[1] == FUNHEAD+4 )
goto NormZero;
3277 else if ( t[1] == *u+FUNHEAD+2 )
goto NormZero;
3280 u = t+FUNHEAD; NEXTARG(u);
3281 if ( *u == -SNUMBER && u[1] == 0 )
goto NormInf;
3283 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3284 else if ( i < ncom-1 && pcom[i+1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3286 if ( AN.PolyNormFlag ) {
3287 if ( AR.PolyFunExp == 0 ) {
3291 else if ( AR.PolyFunExp == 1 ) {
3292 if ( PolyFunMode == 0 ) {
3299 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3305 if ( PolyFunMode == 0 ) {
3312 if ( ExpandRat(BHEAD mmm) != 0 )
3319 if ( AR.PolyFunExp == 0 ) {
3323 else if ( AR.PolyFunExp == 1 ) {
3326 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3333 if ( ExpandRat(BHEAD mmm) != 0 )
3340 else if ( *t > 0 ) {
3341 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3342 && *t != REPLACEMENT && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3347 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3356 if ( ReplaceVeto < 0 ) {
3366 WORD *ma = fillsetexp, *mb, *mc;
3369 if ( *ma != REPLACEMENT ) {
3373 if ( *ma == REPLACEMENT && ReplaceType == -1 ) {
3376 if ( AN.RSsize < 2*ma[1]+SUBEXPSIZE ) {
3377 if ( AN.ReplaceScrat ) M_free(AN.ReplaceScrat,
"AN.ReplaceScrat");
3378 AN.RSsize = 2*ma[1]+SUBEXPSIZE+40;
3379 AN.ReplaceScrat = (WORD *)Malloc1((AN.RSsize+1)*
sizeof(WORD),
"AN.ReplaceScrat");
3382 ReplaceSub = AN.ReplaceScrat;
3383 ReplaceSub += SUBEXPSIZE;
3385 if ( *ma > 0 )
goto NoRep;
3386 if ( *ma <= -FUNCTION ) {
3387 *ReplaceSub++ = FUNTOFUN;
3389 *ReplaceSub++ = -*ma++;
3390 if ( *ma > -FUNCTION )
goto NoRep;
3391 *ReplaceSub++ = -*ma++;
3393 else if ( ma+4 > mb )
goto NoRep;
3395 if ( *ma == -SYMBOL ) {
3396 if ( ma[2] == -SYMBOL && ma+4 <= mb )
3397 *ReplaceSub++ = SYMTOSYM;
3398 else if ( ma[2] == -SNUMBER && ma+4 <= mb ) {
3399 *ReplaceSub++ = SYMTONUM;
3400 if ( ReplaceType == 0 ) {
3401 oldtoprhs = C->numrhs;
3406 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3407 *ReplaceSub++ = SYMTONUM;
3409 *ReplaceSub++ = ma[1];
3418 else if ( ma[2] > 0 ) {
3419 WORD *sstop, *ttstop, n;
3423 while ( ss < sstop ) {
3425 ttstop = tt - ABS(tt[-1]);
3427 while ( ss < ttstop ) {
3428 if ( *ss == INDEX )
goto NoRep;
3434 if ( ReplaceType == 0 ) {
3435 oldtoprhs = C->numrhs;
3439 ss =
AddRHS(AT.ebufnum,1);
3444 while ( --n >= 0 ) *ss++ = *tt++;
3446 C->
rhs[C->numrhs+1] = ss;
3448 *ReplaceSub++ = subtype;
3450 *ReplaceSub++ = ma[1];
3451 *ReplaceSub++ = C->numrhs;
3457 else if ( ( *ma == -VECTOR || *ma == -MINVECTOR ) && ma+4 <= mb ) {
3458 if ( ma[2] == -VECTOR ) {
3459 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOVEC;
3460 else *ReplaceSub++ = VECTOMIN;
3462 else if ( ma[2] == -MINVECTOR ) {
3463 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOMIN;
3464 else *ReplaceSub++ = VECTOVEC;
3470 else if ( ma[2] > 0 ) {
3471 WORD *sstop, *ttstop, *w, *mm, n, count;
3473 if ( *ma == -MINVECTOR ) {
3477 while ( ss < sstop ) {
3486 while ( ss < sstop ) {
3488 ttstop = tt - ABS(tt[-1]);
3491 while ( ss < ttstop ) {
3492 if ( *ss == INDEX ) {
3493 n = ss[1] - 2; ss += 2;
3494 while ( --n >= 0 ) {
3495 if ( *ss < MINSPEC ) count++;
3501 if ( count != 1 )
goto NoRep;
3505 if ( ReplaceType == 0 ) {
3506 oldtoprhs = C->numrhs;
3510 mm =
AddRHS(AT.ebufnum,1);
3511 *ReplaceSub++ = subtype;
3513 *ReplaceSub++ = ma[1];
3514 *ReplaceSub++ = C->numrhs;
3518 while ( (mm + n + 10) > C->
Top )
3520 while ( --n >= 0 ) *mm++ = *w++;
3522 C->
rhs[C->numrhs+1] = mm;
3524 mm =
AddRHS(AT.ebufnum,1);
3528 while ( (mm + n + 13) > C->
Top )
3531 while ( w < sstop ) {
3532 tt = w + *w; ttstop = tt - ABS(tt[-1]);
3534 while ( w < ttstop ) {
3535 if ( *w != INDEX ) {
3544 while ( --n >= 0 ) {
3545 if ( *w >= MINSPEC ) *mm++ = *w++;
3550 if ( n <= 2 ) mm -= 2;
3559 while ( w < tt ) *mm++ = *w++;
3560 *ss = WORDDIF(mm,ss);
3563 C->
rhs[C->numrhs+1] = mm;
3565 if ( mm > C->
Top ) {
3566 MLOCK(ErrorMessageLock);
3567 MesPrint(
"Internal error in Normalize with extra compiler buffer");
3568 MUNLOCK(ErrorMessageLock);
3576 else if ( *ma == -INDEX ) {
3577 if ( ( ma[2] == -INDEX || ma[2] == -VECTOR )
3579 *ReplaceSub++ = INDTOIND;
3580 else if ( ma[1] >= AM.OffsetIndex ) {
3581 if ( ma[2] == -SNUMBER && ma+4 <= mb
3582 && ma[3] >= 0 && ma[3] < AM.OffsetIndex )
3583 *ReplaceSub++ = INDTOIND;
3584 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3585 *ReplaceSub++ = INDTOIND;
3587 *ReplaceSub++ = ma[1];
3598 *ReplaceSub++ = ma[1];
3599 *ReplaceSub++ = ma[3];
3604 AN.ReplaceScrat[1] = ReplaceSub-AN.ReplaceScrat;
3609 while ( mb < m ) *mc++ = *mb++;
3613 if ( ReplaceType > 0 ) {
3614 C->numrhs = oldtoprhs;
3618 if ( ++ReplaceVeto >= 0 )
break;
3629 for ( i = 0; i < neps; i++ ) {
3631 if ( ( t[2] & DIRTYSYMFLAG ) != DIRTYSYMFLAG )
continue;
3632 t[2] &= ~DIRTYSYMFLAG;
3633 if ( AR.Eside == LHSIDE || AR.Eside == LHSIDEX ) {
3642 if ( *r != FUNNYWILD ) { r++;
continue; }
3643 k = r[1]; u = r + 2;
3646 if ( *u != FUNNYWILD ) ncoef = -ncoef;
3649 tt[-2] = FUNNYWILD; tt[-1] = k; m -= 2;
3653 for ( r = t + 1; r < m; r++ ) {
3654 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3655 else if ( *r == *t )
goto NormZero;
3660 for ( r = t + 2; r < tt; r += 2 ) {
3661 if ( r[1] < t[1] ) {
3662 k = r[1]; r[1] = t[1]; t[1] = k; ncoef = -ncoef; }
3663 else if ( r[1] == t[1] )
goto NormZero;
3672 for ( r = t + 1; r < m; r++ ) {
3673 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3674 else if ( *r == *t )
goto NormZero;
3683 for ( i = 0; i < (neps-1); i++ ) {
3685 for ( j = i+1; j < neps; j++ ) {
3687 if ( t[1] > r[1] ) {
3688 peps[i] = m = r; peps[j] = r = t; t = m;
3690 else if ( t[1] == r[1] ) {
3696 m = peps[j]; peps[j] = t; peps[i] = t = m;
3699 else if ( *r++ > *m++ )
break;
3700 }
while ( --k > 0 );
3705 for ( i = 0; i < neps; i++ ) {
3717 for ( i = 0; i < ndel; i += 2, r += 2 ) {
3718 if ( r[1] < r[0] ) { k = *r; *r = r[1]; r[1] = k; }
3720 for ( i = 2; i < ndel; i += 2, t += 2 ) {
3722 for ( j = i; j < ndel; j += 2 ) {
3723 if ( *r > *t ) { r += 2; }
3724 else if ( *r < *t ) {
3725 k = *r; *r++ = *t; *t++ = k;
3726 k = *r; *r++ = *t; *t-- = k;
3729 if ( *++r < t[1] ) {
3730 k = *r; *r = t[1]; t[1] = k;
3748 for ( i = 0; i < nind; i++ ) {
3750 for ( j = i+1; j < nind; j++ ) {
3752 k = *r; *r = *t; *t = k;
3770 for ( i = 2; i < nvec; i += 2 ) {
3772 for ( j = i; j < nvec; j += 2 ) {
3774 if ( *++r < t[1] ) {
3775 k = *r; *r = t[1]; t[1] = k;
3779 else if ( *r < *t ) {
3780 k = *r; *r++ = *t; *t++ = k;
3781 k = *r; *r++ = *t; *t-- = k;
3801 while ( --i >= 0 ) {
3802 if ( *t > t[1] ) { j = *t; *t = t[1]; t[1] = j; }
3808 while ( t < (m-3) ) {
3812 if ( *++r == *++t ) {
3814 if ( ( *r < MAXPOWER && t[1] < MAXPOWER )
3815 || ( *r > -MAXPOWER && t[1] > -MAXPOWER ) ) {
3818 if ( *t > MAXPOWER || *t < -MAXPOWER ) {
3819 MLOCK(ErrorMessageLock);
3820 MesPrint(
"Exponent of dotproduct out of range: %d",*t);
3821 MUNLOCK(ErrorMessageLock);
3837 else if ( *r < *++t ) {
3838 k = *r; *r++ = *t; *t = k;
3843 else if ( *r < *t ) {
3844 k = *r; *r++ = *t; *t++ = k;
3845 k = *r; *r++ = *t; *t = k;
3848 else { r += 2; t--; }
3850 else if ( *r < *t ) {
3851 k = *r; *r++ = *t; *t++ = k;
3852 k = *r; *r++ = *t; *t++ = k;
3853 k = *r; *r++ = *t; *t = k;
3862 if ( ( i = ndot ) > 0 ) {
3877 *m++ = ( i = nsym ) + 2;
3880 if ( t[1] < (2*MAXPOWER) ) {
3881 if ( t[1] & 1 ) { *m++ = 0; *m++ = 1; }
3883 if ( *++t & 2 ) ncoef = -ncoef;
3887 else if ( *t <= NumSymbols && *t > -2*MAXPOWER ) {
3888 if ( ( ( ( t[1] > symbols[*t].maxpower ) && ( symbols[*t].maxpower < MAXPOWER ) ) ||
3889 ( ( t[1] < symbols[*t].minpower ) && ( symbols[*t].minpower > -MAXPOWER ) ) ) &&
3890 ( t[1] < 2*MAXPOWER ) && ( t[1] > -2*MAXPOWER ) ) {
3891 if ( i <= 2 || t[2] != *t )
goto NormZero;
3893 if ( AN.ncmod == 1 && ( AC.modmode & ALSOPOWERS ) != 0 ) {
3894 if ( AC.cmod[0] == 1 ) t[1] = 0;
3895 else if ( t[1] >= 0 ) t[1] = 1 + (t[1]-1)%(AC.cmod[0]-1);
3897 t[1] = -1 - (-t[1]-1)%(AC.cmod[0]-1);
3898 if ( t[1] < 0 ) t[1] += (AC.cmod[0]-1);
3901 if ( ( t[1] < (2*MAXPOWER) && t[1] >= MAXPOWER )
3902 || ( t[1] > -(2*MAXPOWER) && t[1] <= -MAXPOWER ) ) {
3903 MLOCK(ErrorMessageLock);
3904 MesPrint(
"Exponent out of range: %d",t[1]);
3905 MUNLOCK(ErrorMessageLock);
3908 if ( AT.TrimPower && AR.PolyFunVar == *t && t[1] > AR.PolyFunPow ) {
3915 else { *r -= 2; t += 2; }
3918 *m++ = *t++; *m++ = *t++;
3920 }
while ( (i-=2) > 0 ); }
3921 if ( *r <= 2 ) m = r-1;
3927 stop = (WORD *)(((UBYTE *)(termout)) + AM.MaxTer);
3929 if ( ( m + i ) > stop ) {
3930 MLOCK(ErrorMessageLock);
3931 MesPrint(
"Term too complex during normalization");
3932 MUNLOCK(ErrorMessageLock);
3935 if ( ReplaceType >= 0 ) {
3942 if ( ReplaceType == 0 ) {
3943 AT.WorkPointer = termout+*termout;
3944 WildFill(BHEAD term,termout,AN.ReplaceScrat);
3945 termout = term + *term;
3948 AT.WorkPointer = r = termout + *termout;
3949 WildFill(BHEAD r,termout,AN.ReplaceScrat);
3956 r += *term; r -= ABS(r[-1]);
3959 if ( *m >= FUNCTION && m[1] > FUNHEAD &&
3960 functions[*m-FUNCTION].spec != TENSORFUNCTION )
3973 TermFree(n_llnum,
"n_llnum");
3974 TermFree(n_coef,
"NormCoef");
3993 if ( termout < term + *term && termout >= term ) AT.WorkPointer = term + *term;
3994 else AT.WorkPointer = termout;
4000 TermFree(n_llnum,
"n_llnum");
4001 TermFree(n_coef,
"NormCoef");
4005 MLOCK(ErrorMessageLock);
4006 MesPrint(
"Division by zero during normalization");
4007 MUNLOCK(ErrorMessageLock);
4011 MLOCK(ErrorMessageLock);
4012 MesPrint(
"0^0 during normalization of term");
4013 MUNLOCK(ErrorMessageLock);
4017 MLOCK(ErrorMessageLock);
4018 MesPrint(
"0/0 in polyratfun during normalization of term");
4019 MUNLOCK(ErrorMessageLock);
4024 AT.WorkPointer = termout;
4025 TermFree(n_llnum,
"n_llnum");
4026 TermFree(n_coef,
"NormCoef");
4030 TermFree(n_llnum,
"n_llnum");
4031 TermFree(n_coef,
"NormCoef");
4035 MLOCK(ErrorMessageLock);
4037 MUNLOCK(ErrorMessageLock);
4038 TermFree(n_llnum,
"n_llnum");
4039 TermFree(n_coef,
"NormCoef");
4052 WORD ExtraSymbol(WORD sym, WORD pow, WORD nsym, WORD *ppsym, WORD *ncoef)
4060 if ( pow > 2*MAXPOWER || pow < -2*MAXPOWER
4061 || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
4062 MLOCK(ErrorMessageLock);
4063 MesPrint(
"Illegal wildcard power combination.");
4064 MUNLOCK(ErrorMessageLock);
4069 if ( ( sym <= NumSymbols && sym > -MAXPOWER )
4070 && ( symbols[sym].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
4071 *m %= symbols[sym].maxpower;
4072 if ( *m < 0 ) *m += symbols[sym].maxpower;
4073 if ( ( symbols[sym].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
4074 if ( ( ( symbols[sym].maxpower & 1 ) == 0 ) &&
4075 ( *m >= symbols[sym].maxpower/2 ) ) {
4076 *m -= symbols[sym].maxpower/2; *ncoef = -*ncoef;
4081 if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
4082 MLOCK(ErrorMessageLock);
4083 MesPrint(
"Power overflow during normalization");
4084 MUNLOCK(ErrorMessageLock);
4090 { *m = m[2]; m++; *m = m[2]; m++; i++; }
4095 else if ( sym < *m ) {
4103 { m--; m[2] = *m; m--; m[2] = *m; i++; }
4114 WORD DoTheta(PHEAD WORD *t)
4117 WORD k, *r1, *r2, *tstop, type;
4118 WORD ia, *ta, *tb, *stopa, *stopb;
4119 if ( AC.BracketNormalize )
return(-1);
4124 if ( k <= FUNHEAD )
return(1);
4127 if ( r1 == tstop ) {
4131 if ( *t == ARGHEAD ) {
4132 if ( type == THETA )
return(1);
4136 if ( *t == -SNUMBER ) {
4137 if ( t[1] < 0 )
return(0);
4139 if ( type == THETA2 && t[1] == 0 )
return(0);
4146 if ( *t == ABS(k)+1+ARGHEAD ) {
4147 if ( k > 0 )
return(1);
4157 if ( r2 < tstop ) return(-1);
4162 if ( *t == -SNUMBER && *r1 == -SNUMBER ) {
4163 if ( t[1] > r1[1] )
return(0);
4164 else if ( t[1] < r1[1] ) {
4167 else if ( type == THETA )
return(1);
4170 else if ( t[1] == 0 && *t == -SNUMBER ) {
4172 else if ( *t < *r1 )
return(1);
4173 else if ( *t > *r1 )
return(0);
4175 else if ( r1[1] == 0 && *r1 == -SNUMBER ) {
4177 else if ( *t < *r1 )
return(1);
4178 else if ( *t > *r1 )
return(0);
4180 r2 = AT.WorkPointer;
4194 ta += ARGHEAD; tb += ARGHEAD;
4195 while ( ta < stopa ) {
4196 if ( tb >= stopb )
return(0);
4197 if ( ( ia = CompareTerms(ta,tb,(WORD)1) ) < 0 )
return(0);
4198 if ( ia > 0 )
return(1);
4202 if ( type == THETA )
return(1);
4211 WORD DoDelta(WORD *t)
4213 WORD k, *r1, *r2, *tstop, isnum, isnum2, type = *t;
4214 if ( AC.BracketNormalize )
return(-1);
4216 if ( k <= FUNHEAD )
goto argzero;
4217 if ( k == FUNHEAD+ARGHEAD && t[FUNHEAD] == ARGHEAD )
goto argzero;
4224 if ( *t == -SNUMBER ) { isnum = 1; k = t[1]; }
4230 if ( k == *t-ARGHEAD-1 ) isnum = 1;
4234 if ( r1 >= tstop ) {
4235 if ( !isnum )
return(-1);
4236 if ( k == 0 )
goto argzero;
4241 if ( r2 < tstop ) return(-1);
4243 if ( *r1 == -SNUMBER ) { isnum2 = 1; }
4249 if ( k == *r1-ARGHEAD-1 ) isnum2 = 1;
4252 if ( isnum != isnum2 )
return(-1);
4254 while ( t < tstop && r1 < r2 ) {
4256 if ( !isnum )
return(-1);
4261 if ( t != tstop || r1 != r2 ) {
4262 if ( !isnum )
return(-1);
4266 if ( type == DELTA2 )
return(1);
4269 if ( type == DELTA2 )
return(0);
4278 void DoRevert(WORD *fun, WORD *tmp)
4280 WORD *t, *r, *m, *to, *tt, *mm, i, j;
4285 if ( *r == -REVERSEFUNCTION ) {
4287 while ( mm < to ) *m++ = *mm++;
4290 fun[2] |= DIRTYSYMFLAG;
4292 else if ( *r <= -FUNCTION ) r++;
4294 if ( *r == -INDEX && r[1] < MINSPEC ) *r = -VECTOR;
4299 if ( ( *r > ARGHEAD )
4300 && ( r[ARGHEAD+1] == REVERSEFUNCTION )
4301 && ( *r == (r[ARGHEAD]+ARGHEAD) )
4302 && ( r[ARGHEAD] == (r[ARGHEAD+2]+4) )
4303 && ( *(r+*r-3) == 1 )
4304 && ( *(r+*r-2) == 1 )
4305 && ( *(r+*r-1) == 3 ) ) {
4317 while ( --j >= 0 ) {
4320 while ( --i >= 0 ) {
4327 else if ( *t <= -FUNCTION ) *m++ = *t++;
4328 else { *m++ = *t++; *m++ = *t++; }
4339 fun[1] = WORDDIF(t,fun);
4341 fun[2] |= DIRTYSYMFLAG;
4364 #define MAXNUMBEROFNONCOMTERMS 2
4366 WORD DetCommu(WORD *terms)
4368 WORD *t, *tnext, *tstop;
4370 if ( *terms == 0 )
return(0);
4371 if ( terms[*terms] == 0 )
return(0);
4375 tstop = tnext - ABS(tnext[-1]);
4377 while ( t < tstop ) {
4378 if ( *t >= FUNCTION ) {
4379 if ( functions[*t-FUNCTION].commute ) {
4381 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4385 else if ( *t == SUBEXPRESSION ) {
4386 if ( cbuf[t[4]].CanCommu[t[2]] ) {
4388 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4392 else if ( *t == EXPRESSION ) {
4394 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4397 else if ( *t == DOLLAREXPRESSION ) {
4404 if ( cbuf[AM.dbufnum].CanCommu[t[2]] ) {
4406 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4425 WORD DoesCommu(WORD *term)
4429 if ( *term == 0 )
return(0);
4430 tstop = term + *term;
4431 tstop = tstop - ABS(tstop[-1]);
4433 while ( term < tstop ) {
4434 if ( ( *term >= FUNCTION ) && ( functions[*term-FUNCTION].commute ) ) {
4436 if ( num >= MAXNUMBEROFNONCOMTERMS )
return(num);
4451 WORD *PolyNormPoly (PHEAD WORD *Poly) {
4454 WORD *buffer = AT.WorkPointer;
4456 if (
NewSort(BHEAD0) ) { Terminate(-1); }
4468 if (
EndSort(BHEAD buffer,1) < 0 ) {
4473 while ( *p ) p += *p;
4475 AT.WorkPointer = p + 1;
4502 WORD *EvaluateGcd(PHEAD WORD *subterm)
4505 WORD *oldworkpointer = AT.WorkPointer, *work1, *work2, *work3;
4506 WORD *t, *tt, *ttt, *t1, *t2, *t3, *t4, *tstop;
4509 WORD *lnum=n_llnum+1;
4510 WORD *num1, *num2, *num3, *den1, *den2, *den3;
4511 WORD sizenum1, sizenum2, sizenum3, sizeden1, sizeden2, sizeden3;
4512 int i, isnumeric = 0, numarg = 0 ;
4518 tt = subterm + subterm[1]; t = subterm + FUNHEAD;
4522 if ( *t == -SNUMBER ) {
4525 MLOCK(ErrorMessageLock);
4526 MesPrint(
"Trying to take the GCD involving a zero term.");
4527 MUNLOCK(ErrorMessageLock);
4531 t1 = subterm + FUNHEAD;
4532 while ( gcdnum > 1 && t1 < tt ) {
4533 if ( *t1 == -SNUMBER ) {
4535 if ( stor == 0 )
goto gcdzero;
4536 if ( GcdLong(BHEAD (UWORD *)&stor,1,(UWORD *)&gcdnum,1,
4537 (UWORD *)lnum,&nnum) )
goto FromGCD;
4542 else if ( *t1 == -SYMBOL )
goto gcdisone;
4543 else if ( *t1 < 0 )
goto gcdillegal;
4549 ct = *ttt; *ttt = 0;
4551 t1 = PolyNormPoly(BHEAD t1+ARGHEAD);
4560 while ( t3 > t2 && *t3 == 0 ) { t3--; i--; }
4561 if ( GcdLong(BHEAD (UWORD *)t2,(WORD)i,(UWORD *)&gcdnum,1,
4562 (UWORD *)lnum,&nnum) ) {
4567 if ( gcdnum == 1 ) {
4574 AT.WorkPointer = oldworkpointer;
4576 if ( gcdnum == 1 )
goto gcdisone;
4577 oldworkpointer[0] = 4;
4578 oldworkpointer[1] = gcdnum;
4579 oldworkpointer[2] = 1;
4580 oldworkpointer[3] = 3;
4581 oldworkpointer[4] = 0;
4582 AT.WorkPointer = oldworkpointer + 5;
4583 return(oldworkpointer);
4585 else if ( *t == -SYMBOL ) {
4586 t1 = subterm + FUNHEAD;
4589 if ( *t1 == -SNUMBER )
goto gcdisone;
4590 if ( *t1 == -SYMBOL ) {
4591 if ( t1[1] != i )
goto gcdisone;
4594 if ( *t1 < 0 )
goto gcdillegal;
4596 ct = *ttt; *ttt = 0;
4598 t2 = PolyNormPoly(BHEAD t1+ARGHEAD);
4600 else t2 = t1 + ARGHEAD;
4604 tstop = t2 - ABS(t2[-1]);
4605 while ( t3 < tstop ) {
4606 if ( *t3 != SYMBOL ) {
4613 if ( *t4 == i && t4[1] > 0 )
goto nextterminarg;
4623 AT.WorkPointer = oldworkpointer;
4625 oldworkpointer[0] = 8;
4626 oldworkpointer[1] = SYMBOL;
4627 oldworkpointer[2] = 4;
4628 oldworkpointer[3] = t[1];
4629 oldworkpointer[4] = 1;
4630 oldworkpointer[5] = 1;
4631 oldworkpointer[6] = 1;
4632 oldworkpointer[7] = 3;
4633 oldworkpointer[8] = 0;
4634 AT.WorkPointer = oldworkpointer+9;
4635 return(oldworkpointer);
4637 else if ( *t < 0 ) {
4639 MLOCK(ErrorMessageLock);
4640 MesPrint(
"Illegal object in gcd_ function. Object not a number or a symbol.");
4641 MUNLOCK(ErrorMessageLock);
4644 else if ( ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4645 else if ( t[1] != 0 ) {
4646 ttt = t + *t; ct = *ttt; *ttt = 0;
4647 t = PolyNormPoly(BHEAD t+ARGHEAD);
4649 if ( t[*t] == 0 && ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4650 AT.WorkPointer = oldworkpointer;
4665 AT.WorkPointer = oldworkpointer;
4667 t = subterm + FUNHEAD;
4668 for ( i = 1; i < isnumeric; i++ ) {
4672 ttt = t + *t; ct = *ttt; *ttt = 0;
4673 t = PolyNormPoly(BHEAD t+ARGHEAD);
4677 i = (ABS(t[-1])-1)/2;
4680 sizenum1 = sizeden1 = i;
4681 while ( sizenum1 > 1 && num1[sizenum1-1] == 0 ) sizenum1--;
4682 while ( sizeden1 > 1 && den1[sizeden1-1] == 0 ) sizeden1--;
4683 work1 = AT.WorkPointer+1; work2 = work1+sizenum1;
4684 for ( i = 0; i < sizenum1; i++ ) work1[i] = num1[i];
4685 for ( i = 0; i < sizeden1; i++ ) work2[i] = den1[i];
4686 num1 = work1; den1 = work2;
4687 AT.WorkPointer = work2 = work2 + sizeden1;
4688 t = subterm + FUNHEAD;
4690 ttt = t + *t; ct = *ttt; *ttt = 0;
4692 t = PolyNormPoly(BHEAD t+ARGHEAD);
4697 i = (ABS(t[-1])-1)/2;
4700 sizenum2 = sizeden2 = i;
4701 while ( sizenum2 > 1 && num2[sizenum2-1] == 0 ) sizenum2--;
4702 while ( sizeden2 > 1 && den2[sizeden2-1] == 0 ) sizeden2--;
4703 num3 = AT.WorkPointer;
4704 if ( GcdLong(BHEAD (UWORD *)num2,sizenum2,(UWORD *)num1,sizenum1,
4705 (UWORD *)num3,&sizenum3) )
goto FromGCD;
4706 sizenum1 = sizenum3;
4707 for ( i = 0; i < sizenum1; i++ ) num1[i] = num3[i];
4708 den3 = AT.WorkPointer;
4709 if ( GcdLong(BHEAD (UWORD *)den2,sizeden2,(UWORD *)den1,sizeden1,
4710 (UWORD *)den3,&sizeden3) )
goto FromGCD;
4711 sizeden1 = sizeden3;
4712 for ( i = 0; i < sizeden1; i++ ) den1[i] = den3[i];
4713 if ( sizenum1 == 1 && num1[0] == 1 && sizeden1 == 1 && den1[1] == 1 )
4718 AT.WorkPointer = work2;
4720 AT.WorkPointer = oldworkpointer;
4724 if ( sizenum1 > sizeden1 ) {
4725 while ( sizenum1 > sizeden1 ) den1[sizeden1++] = 0;
4727 else if ( sizenum1 < sizeden1 ) {
4728 while ( sizenum1 < sizeden1 ) num1[sizenum1++] = 0;
4733 if ( num1 != t ) { NCOPY(t,num1,sizenum1); }
4735 if ( den1 != t ) { NCOPY(t,den1,sizeden1); }
4740 return(oldworkpointer);
4747 t = subterm + FUNHEAD;
4748 AT.WorkPointer += AM.MaxTer/
sizeof(WORD);
4749 work2 = AT.WorkPointer;
4756 work1 = AT.WorkPointer;
4757 ttt = t + *t; ct = *ttt; *ttt = 0;
4758 t = PolyNormPoly(BHEAD t+ARGHEAD);
4759 if ( *work1 < AT.WorkPointer-work1 ) {
4768 *AT.WorkPointer++ = 0;
4774 if ( work2 != work3 ) {
4775 work1 = PolyGCD2(BHEAD work1,work2);
4777 while ( *work2 ) work2 += *work2;
4781 while ( *work2 ) work2 += *work2;
4782 size = work2 - work1 + 1;
4784 NCOPY(t,work1,size);
4786 return(oldworkpointer);
4789 oldworkpointer[0] = 4;
4790 oldworkpointer[1] = 1;
4791 oldworkpointer[2] = 1;
4792 oldworkpointer[3] = 3;
4793 oldworkpointer[4] = 0;
4794 AT.WorkPointer = oldworkpointer+5;
4795 return(oldworkpointer);
4797 MLOCK(ErrorMessageLock);
4798 MesCall(
"EvaluateGcd");
4799 MUNLOCK(ErrorMessageLock);
4816 int TreatPolyRatFun(PHEAD WORD *prf)
4818 WORD *t, *tstop, *r, *rstop, *m, *mstop;
4819 WORD exp1 = MAXPOWER, exp2 = MAXPOWER;
4822 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
4823 if ( exp1 > 1 ) exp1 = 1;
4827 if ( exp1 > 0 ) exp1 = 0;
4834 while ( t < tstop ) {
4840 rstop = t - ABS(t[-1]);
4841 while ( r < rstop ) {
4842 if ( *r != SYMBOL ) { r += r[1];
continue; }
4846 while ( m < mstop ) {
4847 if ( *m == AR.PolyFunVar ) {
4848 if ( m[1] < exp1 ) exp1 = m[1];
4854 if ( exp1 > 0 ) exp1 = 0;
4859 if ( exp1 > 0 ) exp1 = 0;
4865 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
4866 if ( exp2 > 1 ) exp2 = 1;
4869 if ( exp2 > 0 ) exp2 = 0;
4875 while ( t < tstop ) {
4881 rstop = t - ABS(t[-1]);
4882 while ( r < rstop ) {
4883 if ( *r != SYMBOL ) { r += r[1];
continue; }
4887 while ( m < mstop ) {
4888 if ( *m == AR.PolyFunVar ) {
4889 if ( m[1] < exp2 ) exp2 = m[1];
4895 if ( exp2 > 0 ) exp2 = 0;
4900 if ( exp2 > 0 ) exp2 = 0;
4913 *t++ = -SNUMBER; *t++ = 1;
4914 *t++ = -SNUMBER; *t++ = 1;
4916 else if ( exp1 > 0 ) {
4918 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
4924 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
4925 *t++ = exp1; *t++ = 1; *t++ = 1; *t++ = 3;
4927 *t++ = -SNUMBER; *t++ = 1;
4930 *t++ = -SNUMBER; *t++ = 1;
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;
4952 void DropCoefficient(PHEAD WORD *term)
4955 WORD *t = term + *term;
4957 n = t[-1]; na = ABS(n);
4959 if ( n == 3 && t[0] == 1 && t[1] == 1 )
return;
4961 t[0] = 1; t[1] = 1; t[2] = 3;
4970 void DropSymbols(PHEAD WORD *term)
4973 WORD *tend = term + *term, *t1, *t2, *tstop;
4974 tstop = tend - ABS(tend[-1]);
4976 while ( t1 < tstop ) {
4977 if ( *t1 == SYMBOL ) {
4980 while ( t2 < tend ) *t1++ = *t2++;
5003 WORD buffer[7*NORMSIZE], *t, *b, *bb, *tt, *m, *tstop;
5006 *b++ = SYMBOL; *b++ = 2;
5008 tstop = t - ABS(t[-1]);
5010 while ( t < tstop ) {
5011 if ( *t == SYMBOL && t < tstop ) {
5012 for ( i = 2; i < t[1]; i += 2 ) {
5015 if ( bb[0] == t[i] ) {
5017 if ( bb[1] > MAXPOWER || bb[1] < -MAXPOWER ) {
5018 MLOCK(ErrorMessageLock);
5019 MesPrint(
"Power in SymbolNormalize out of range");
5020 MUNLOCK(ErrorMessageLock);
5026 bb[0] = bb[2]; bb[1] = bb[3]; bb += 2;
5031 else if ( bb[0] > t[i] ) {
5033 while ( m > bb ) { m[1] = m[-1]; m[0] = m[-2]; m -= 2; }
5042 *b++ = t[i]; *b++ = t[i+1];
5048 MLOCK(ErrorMessageLock);
5049 MesPrint(
"Illegal term in SymbolNormalize");
5050 MUNLOCK(ErrorMessageLock);
5055 buffer[1] = b - buffer;
5059 if ( AT.LeaveNegative == 0 ) {
5060 b = buffer; bb = b + b[1]; b += 3;
5063 MLOCK(ErrorMessageLock);
5064 MesPrint(
"Negative power in SymbolNormalize");
5065 MUNLOCK(ErrorMessageLock);
5077 b = buffer; tt = term + 1;
5078 if ( i > 2 ) { NCOPY(tt,b,i) }
5081 if ( i < 0 ) i = -i;
5082 *term -= (tstop-tt);
5096 int TestFunFlag(PHEAD WORD *tfun)
5098 WORD *t, *tstop, *r, *rstop, *m, *mstop;
5099 if ( functions[*tfun-FUNCTION].spec )
return(0);
5100 tstop = tfun + tfun[1];
5102 while ( t < tstop ) {
5103 if ( *t < 0 ) { NEXTARG(t);
continue; }
5105 if ( t[1] == 0 ) { t = rstop;
continue; }
5107 while ( r < rstop ) {
5108 m = r+1; mstop = r+*r; mstop -= ABS(mstop[-1]);
5109 while ( m < mstop ) {
5110 if ( *m == SUBEXPRESSION || *m == EXPRESSION || *m == DOLLAREXPRESSION )
return(1);
5111 if ( ( *m >= FUNCTION ) && ( ( m[2] & DIRTYFLAG ) == DIRTYFLAG )
5112 && ( *m != REPLACEMENT ) && TestFunFlag(BHEAD m) )
return(1);
5127 WORD BracketNormalize(PHEAD WORD *term)
5129 WORD *stop = term+*term-3, *t, *tt, *tstart, *r;
5130 WORD *oldwork = AT.WorkPointer;
5133 termout = AT.WorkPointer = term+*term;
5137 tt = termout+1; t = term+1;
5138 while ( t < stop ) {
5139 if ( *t >= FUNCTION ) { i = t[1]; NCOPY(tt,t,i); }
5142 if ( tt > termout+1 && tt-termout-1 > termout[2] ) {
5143 r = termout+1; ii = tt-r;
5144 for ( i = 0; i < ii-FUNHEAD; i += FUNHEAD ) {
5145 for ( j = i+FUNHEAD; j > 0; j -= FUNHEAD ) {
5146 if ( functions[r[j-FUNHEAD]-FUNCTION].commute
5147 && functions[r[j]-FUNCTION].commute == 0 )
break;
5148 if ( r[j-FUNHEAD] > r[j] ) EXCH(r[j-FUNHEAD],r[j])
5154 tstart = tt; t = term + 1; *tt++ = DELTA; *tt++ = 2;
5155 while ( t < stop ) {
5156 if ( *t == DELTA ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5159 if ( tstart[1] > 2 ) {
5160 for ( r = tstart+2; r < tstart+tstart[1]; r += 2 ) {
5161 if ( r[0] > r[1] ) EXCH(r[0],r[1])
5164 if ( tstart[1] > 4 ) {
5165 r = tstart+2; ii = tstart[1]-2;
5166 for ( i = 0; i < ii-2; i += 2 ) {
5167 for ( j = i+2; j > 0; j -= 2 ) {
5168 if ( r[j-2] > r[j] ) {
5172 else if ( r[j-2] < r[j] )
break;
5174 if ( r[j-1] > r[j+1] ) EXCH(r[j-1],r[j+1])
5179 tt = tstart+tstart[1];
5181 else if ( tstart[1] == 2 ) { tt = tstart; }
5184 tstart = tt; t = term + 1; *tt++ = INDEX; *tt++ = 2;
5185 while ( t < stop ) {
5186 if ( *t == INDEX ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5189 if ( tstart[1] >= 4 ) {
5190 r = tstart+2; ii = tstart[1]-2;
5191 for ( i = 0; i < ii-1; i += 1 ) {
5192 for ( j = i+1; j > 0; j -= 1 ) {
5193 if ( r[j-1] > r[j] ) EXCH(r[j-1],r[j])
5197 tt = tstart+tstart[1];
5199 else if ( tstart[1] == 2 ) { tt = tstart; }
5202 tstart = tt; t = term + 1; *tt++ = DOTPRODUCT; *tt++ = 2;
5203 while ( t < stop ) {
5204 if ( *t == DOTPRODUCT ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5207 if ( tstart[1] > 5 ) {
5208 r = tstart+2; ii = tstart[1]-2;
5209 for ( i = 0; i < ii; i += 3 ) {
5210 if ( r[i] > r[i+1] ) EXCH(r[i],r[i+1])
5212 for ( i = 0; i < ii-3; i += 3 ) {
5213 for ( j = i+3; j > 0; j -= 3 ) {
5214 if ( r[j-3] < r[j] )
break;
5215 if ( r[j-3] > r[j] ) {
5220 if ( r[j-2] > r[j+1] ) EXCH(r[j-2],r[j+1])
5225 tt = tstart+tstart[1];
5227 else if ( tstart[1] == 2 ) { tt = tstart; }
5229 if ( tstart[2] > tstart[3] ) EXCH(tstart[2],tstart[3])
5233 tstart = tt; t = term + 1; *tt++ = SYMBOL; *tt++ = 2;
5234 while ( t < stop ) {
5235 if ( *t == SYMBOL ) { 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] ) EXCH(r[j-2],r[j])
5246 tt = tstart+tstart[1];
5248 else if ( tstart[1] == 2 ) { tt = tstart; }
5251 tstart = tt; t = term + 1; *tt++ = SETSET; *tt++ = 2;
5252 while ( t < stop ) {
5253 if ( *t == SETSET ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5256 if ( tstart[1] > 4 ) {
5257 r = tstart+2; ii = tstart[1]-2;
5258 for ( i = 0; i < ii-2; i += 2 ) {
5259 for ( j = i+2; j > 0; j -= 2 ) {
5260 if ( r[j-2] > r[j] ) {
5267 tt = tstart+tstart[1];
5269 else if ( tstart[1] == 2 ) { tt = tstart; }
5271 *tt++ = 1; *tt++ = 1; *tt++ = 3;
5272 t = term; i = *termout = tt - termout; tt = termout;
5274 AT.WorkPointer = oldwork;
WORD Compare1(WORD *, WORD *, WORD)
WORD CompareSymbols(WORD *, WORD *, WORD)
int SymbolNormalize(WORD *)
WORD StoreTerm(PHEAD WORD *)
WORD NextPrime(PHEAD WORD)
WORD CompCoef(WORD *, WORD *)
LONG EndSort(PHEAD WORD *, int)