FORM  4.1
normal.c
Go to the documentation of this file.
1 
10 /* #[ License : */
11 /*
12  * Copyright (C) 1984-2013 J.A.M. Vermaseren
13  * When using this file you are requested to refer to the publication
14  * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
15  * This is considered a matter of courtesy as the development was paid
16  * for by FOM the Dutch physics granting agency and we would like to
17  * be able to track its scientific use to convince FOM of its value
18  * for the community.
19  *
20  * This file is part of FORM.
21  *
22  * FORM is free software: you can redistribute it and/or modify it under the
23  * terms of the GNU General Public License as published by the Free Software
24  * Foundation, either version 3 of the License, or (at your option) any later
25  * version.
26  *
27  * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
28  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
29  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
30  * details.
31  *
32  * You should have received a copy of the GNU General Public License along
33  * with FORM. If not, see <http://www.gnu.org/licenses/>.
34  */
35 /* #] License : */
36 /*
37  #[ Includes : normal.c
38 */
39 
40 #include "form3.h"
41 
42 /*
43  #] Includes :
44  #[ Normalize :
45  #[ Commute :
46 
47  This function gets two adjacent function pointers and decides
48  whether these two functions should be exchanged to obtain a
49  natural ordering.
50 
51  Currently there is only an ordering of gamma matrices belonging
52  to different spin lines.
53 
54 */
55 
56 WORD Commute(WORD *fleft, WORD *fright)
57 {
58  if ( *fleft >= GAMMA && *fleft <= GAMMASEVEN
59  && *fright >= GAMMA && *fright <= GAMMASEVEN ) {
60  if ( fleft[FUNHEAD] < AM.OffsetIndex && fleft[FUNHEAD] > fright[FUNHEAD] )
61  return(1);
62  }
63 /*
64  if other conditions will come here, keep in mind that if *fleft < 0
65  or *fright < 0 they are arguments in the exponent function as in f^2
66 */
67  return(0);
68 }
69 
70 /*
71  #] Commute :
72  #[ Normalize :
73 
74  This is the big normalization routine. It has a great need
75  to be economical.
76  There is a fixed limit to the number of objects coming in.
77  Something should be done about it.
78 
79 */
80 
81 WORD Normalize(PHEAD WORD *term)
82 {
83 /*
84  #[ Declarations :
85 */
86  GETBIDENTITY
87  WORD *t, *m, *r, i, j, k, l, nsym, *ss, *tt, *u;
88  WORD shortnum, stype;
89  WORD *stop, *to = 0, *from = 0;
90 /*
91  The next variables would be better off in the AT.WorkSpace (?)
92  or as static global variables. Now they make stackallocations
93  rather bothersome.
94 */
95  WORD psym[7*NORMSIZE],*ppsym;
96  WORD pvec[NORMSIZE],*ppvec,nvec;
97  WORD pdot[3*NORMSIZE],*ppdot,ndot;
98  WORD pdel[2*NORMSIZE],*ppdel,ndel;
99  WORD pind[NORMSIZE],nind;
100  WORD *peps[NORMSIZE/3],neps;
101  WORD *pden[NORMSIZE/3],nden;
102  WORD *pcom[NORMSIZE],ncom;
103  WORD *pnco[NORMSIZE],nnco;
104  WORD *pcon[2*NORMSIZE],ncon; /* Pointer to contractable indices */
105  WORD ncoef; /* Accumulator for the coefficient */
106  WORD *lnum=AT.n_llnum+1,nnum; /* Scratch for factorials */
107  WORD *termout, oldtoprhs = 0, subtype;
108  WORD ReplaceType, didcontr, regval = 0;
109  WORD *ReplaceSub;
110  WORD *fillsetexp;
111  CBUF *C = cbuf+AT.ebufnum;
112  WORD *ANsc = 0, *ANsm = 0, *ANsr = 0;
113  LONG oldcpointer = 0;
114 /*
115  int termflag;
116 */
117 /*
118  #] Declarations :
119  #[ Setup :
120 PrintTerm(term,"Normalize");
121 */
122 
123 Restart:
124  didcontr = 0;
125  ReplaceType = -1;
126  t = term;
127  if ( !*t ) return(regval);
128  r = t + *t;
129  ncoef = r[-1];
130  i = ABS(ncoef);
131  r -= i;
132  m = r;
133  t = AT.n_coef;
134  NCOPY(t,r,i);
135  termout = AT.WorkPointer;
136  AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
137  fillsetexp = termout+1;
138  AN.PolyNormFlag = 0;
139 /*
140  termflag = 0;
141 */
142 /*
143  #] Setup :
144  #[ First scan :
145 */
146  nsym = nvec = ndot = ndel = neps = nden =
147  nind = ncom = nnco = ncon = 0;
148  ppsym = psym;
149  ppvec = pvec;
150  ppdot = pdot;
151  ppdel = pdel;
152  t = term + 1;
153 conscan:;
154  if ( t < m ) do {
155  r = t + t[1];
156  switch ( *t ) {
157  case SYMBOL :
158  t += 2;
159  from = m;
160  do {
161  if ( t[1] == 0 ) {
162 /* if ( *t == 0 || *t == MAXPOWER ) goto NormZZ; */
163  t += 2;
164  goto NextSymbol;
165  }
166  if ( *t <= DENOMINATORSYMBOL && *t >= COEFFSYMBOL ) {
167 /*
168  if ( AN.NoScrat2 == 0 ) {
169  AN.NoScrat2 = (UWORD *)Malloc1((AM.MaxTal+2)*sizeof(UWORD),"Normalize");
170  }
171 */
172  if ( AN.cTerm ) m = AN.cTerm;
173  else m = term;
174  m += *m;
175  ncoef = REDLENG(ncoef);
176  if ( *t == COEFFSYMBOL ) {
177  i = t[1];
178  nnum = REDLENG(m[-1]);
179  m -= ABS(m[-1]);
180  if ( i > 0 ) {
181  while ( i > 0 ) {
182  if ( MulRat(BHEAD (UWORD *)AT.n_coef,ncoef,(UWORD *)m,nnum,
183  (UWORD *)AT.n_coef,&ncoef) ) goto FromNorm;
184  i--;
185  }
186  }
187  else if ( i < 0 ) {
188  while ( i < 0 ) {
189  if ( DivRat(BHEAD (UWORD *)AT.n_coef,ncoef,(UWORD *)m,nnum,
190  (UWORD *)AT.n_coef,&ncoef) ) goto FromNorm;
191  i++;
192  }
193  }
194  }
195  else {
196  i = m[-1];
197  nnum = (ABS(i)-1)/2;
198  if ( *t == NUMERATORSYMBOL ) { m -= nnum + 1; }
199  else { m--; }
200  while ( *m == 0 && nnum > 1 ) { m--; nnum--; }
201  m -= nnum;
202  if ( i < 0 && *t == NUMERATORSYMBOL ) nnum = -nnum;
203  i = t[1];
204  if ( i > 0 ) {
205  while ( i > 0 ) {
206  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)m,nnum) )
207  goto FromNorm;
208  i--;
209  }
210  }
211  else if ( i < 0 ) {
212  while ( i < 0 ) {
213  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)m,nnum) )
214  goto FromNorm;
215  i++;
216  }
217  }
218  }
219  ncoef = INCLENG(ncoef);
220  t += 2;
221  goto NextSymbol;
222  }
223  else if ( *t == DIMENSIONSYMBOL ) {
224  if ( AN.cTerm ) m = AN.cTerm;
225  else m = term;
226  k = DimensionTerm(m);
227  if ( k == 0 ) goto NormZero;
228  if ( k == MAXPOSITIVE ) {
229  MLOCK(ErrorMessageLock);
230  MesPrint("Dimension_ is undefined in term %t");
231  MUNLOCK(ErrorMessageLock);
232  goto NormMin;
233  }
234  if ( k == -MAXPOSITIVE ) {
235  MLOCK(ErrorMessageLock);
236  MesPrint("Dimension_ out of range in term %t");
237  MUNLOCK(ErrorMessageLock);
238  goto NormMin;
239  }
240  if ( k > 0 ) { *((UWORD *)lnum) = k; nnum = 1; }
241  else { *((UWORD *)lnum) = -k; nnum = -1; }
242  ncoef = REDLENG(ncoef);
243  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) ) goto FromNorm;
244  ncoef = INCLENG(ncoef);
245  t += 2;
246  goto NextSymbol;
247  }
248  if ( ( *t >= MAXPOWER && *t < 2*MAXPOWER )
249  || ( *t < -MAXPOWER && *t > -2*MAXPOWER ) ) {
250 /*
251  #[ TO SNUMBER :
252 */
253  if ( *t < 0 ) {
254  *t += MAXPOWER;
255  *t = -*t;
256  if ( t[1] & 1 ) ncoef = -ncoef;
257  }
258  else if ( *t == MAXPOWER ) {
259  if ( t[1] > 0 ) goto NormZero;
260  goto NormInf;
261  }
262  else {
263  *t -= MAXPOWER;
264  }
265  lnum[0] = *t;
266  nnum = 1;
267  if ( t[1] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[1]))) )
268  goto FromNorm;
269  ncoef = REDLENG(ncoef);
270  if ( t[1] < 0 ) {
271  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) )
272  goto FromNorm;
273  }
274  else if ( t[1] > 0 ) {
275  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) )
276  goto FromNorm;
277  }
278  ncoef = INCLENG(ncoef);
279 /*
280  #] TO SNUMBER :
281 */
282  t += 2;
283  goto NextSymbol;
284  }
285  if ( ( *t <= NumSymbols && *t > -MAXPOWER )
286  && ( symbols[*t].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
287  if ( t[1] <= 2*MAXPOWER && t[1] >= -2*MAXPOWER ) {
288  t[1] %= symbols[*t].maxpower;
289  if ( t[1] < 0 ) t[1] += symbols[*t].maxpower;
290  if ( ( symbols[*t].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
291  if ( ( ( symbols[*t].maxpower & 1 ) == 0 ) &&
292  ( t[1] >= symbols[*t].maxpower/2 ) ) {
293  t[1] -= symbols[*t].maxpower/2; ncoef = -ncoef;
294  }
295  }
296  if ( t[1] == 0 ) { t += 2; goto NextSymbol; }
297  }
298  }
299  i = nsym;
300  m = ppsym;
301  if ( i > 0 ) do {
302  m -= 2;
303  if ( *t == *m ) {
304  t++; m++;
305  if ( *t > 2*MAXPOWER || *t < -2*MAXPOWER
306  || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
307  MLOCK(ErrorMessageLock);
308  MesPrint("Illegal wildcard power combination.");
309  MUNLOCK(ErrorMessageLock);
310  goto NormMin;
311  }
312  *m += *t;
313  if ( ( t[-1] <= NumSymbols && t[-1] > -MAXPOWER )
314  && ( symbols[t[-1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
315  *m %= symbols[t[-1]].maxpower;
316  if ( *m < 0 ) *m += symbols[t[-1]].maxpower;
317  if ( ( symbols[t[-1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
318  if ( ( ( symbols[t[-1]].maxpower & 1 ) == 0 ) &&
319  ( *m >= symbols[t[-1]].maxpower/2 ) ) {
320  *m -= symbols[t[-1]].maxpower/2; ncoef = -ncoef;
321  }
322  }
323  }
324  if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
325  MLOCK(ErrorMessageLock);
326  MesPrint("Power overflow during normalization");
327  MUNLOCK(ErrorMessageLock);
328  goto NormMin;
329  }
330  if ( !*m ) {
331  m--;
332  while ( i < nsym )
333  { *m = m[2]; m++; *m = m[2]; m++; i++; }
334  ppsym -= 2;
335  nsym--;
336  }
337  t++;
338  goto NextSymbol;
339  }
340  } while ( *t < *m && --i > 0 );
341  m = ppsym;
342  while ( i < nsym )
343  { m--; m[2] = *m; m--; m[2] = *m; i++; }
344  *m++ = *t++;
345  *m = *t++;
346  ppsym += 2;
347  nsym++;
348 NextSymbol:;
349  } while ( t < r );
350  m = from;
351  break;
352  case VECTOR :
353  t += 2;
354  do {
355  if ( t[1] == FUNNYVEC ) {
356  pind[nind++] = *t;
357  t += 2;
358  }
359  else if ( t[1] < 0 ) {
360  if ( *t == NOINDEX && t[1] == NOINDEX ) t += 2;
361  else {
362  *ppdot++ = *t++; *ppdot++ = *t++; *ppdot++ = 1; ndot++;
363  }
364  }
365  else { *ppvec++ = *t++; *ppvec++ = *t++; nvec += 2; }
366  } while ( t < r );
367  break;
368  case DOTPRODUCT :
369  t += 2;
370  do {
371  if ( t[2] == 0 ) t += 3;
372  else if ( ndot > 0 && t[0] == ppdot[-3]
373  && t[1] == ppdot[-2] ) {
374  ppdot[-1] += t[2];
375  t += 3;
376  if ( ppdot[-1] == 0 ) { ppdot -= 3; ndot--; }
377  }
378  else {
379  *ppdot++ = *t++; *ppdot++ = *t++;
380  *ppdot++ = *t++; ndot++;
381  }
382  } while ( t < r );
383  break;
384  case HAAKJE :
385  break;
386  case SETSET:
387  if ( WildFill(BHEAD termout,term,AT.dummysubexp) < 0 ) goto FromNorm;
388  i = *termout;
389  t = termout; m = term;
390  NCOPY(m,t,i);
391  goto Restart;
392  case DOLLAREXPRESSION :
393 /*
394  We have DOLLAREXPRESSION,4,number,power
395  Replace by SUBEXPRESSION and exit elegantly to let
396  TestSub pick it up. Of course look for special cases first.
397  Note that we have a special compiler buffer for the values.
398 */
399  if ( AR.Eside != LHSIDE ) {
400  DOLLARS d = Dollars + t[2];
401 #ifdef WITHPTHREADS
402  int nummodopt, ptype = -1;
403  if ( AS.MultiThreaded ) {
404  for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
405  if ( t[2] == ModOptdollars[nummodopt].number ) break;
406  }
407  if ( nummodopt < NumModOptdollars ) {
408  ptype = ModOptdollars[nummodopt].type;
409  if ( ptype == MODLOCAL ) {
410  d = ModOptdollars[nummodopt].dstruct+AT.identity;
411  }
412  else {
413  LOCK(d->pthreadslockread);
414  }
415  }
416  }
417 #endif
418  if ( d->type == DOLZERO ) {
419 #ifdef WITHPTHREADS
420  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
421 #endif
422  if ( t[3] == 0 ) goto NormZZ;
423  if ( t[3] < 0 ) goto NormInf;
424  goto NormZero;
425  }
426  else if ( d->type == DOLNUMBER ) {
427  nnum = d->where[0];
428  if ( nnum > 0 ) {
429  nnum = d->where[nnum-1];
430  if ( nnum < 0 ) { ncoef = -ncoef; nnum = -nnum; }
431  nnum = (nnum-1)/2;
432  for ( i = 1; i <= nnum; i++ ) lnum[i-1] = d->where[i];
433  }
434  if ( nnum == 0 || ( nnum == 1 && lnum[0] == 0 ) ) {
435 #ifdef WITHPTHREADS
436  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
437 #endif
438  if ( t[3] < 0 ) goto NormInf;
439  else if ( t[3] == 0 ) goto NormZZ;
440  goto NormZero;
441  }
442  if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) ) goto FromNorm;
443  ncoef = REDLENG(ncoef);
444  if ( t[3] < 0 ) {
445  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) ) {
446 #ifdef WITHPTHREADS
447  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
448 #endif
449  goto FromNorm;
450  }
451  }
452  else if ( t[3] > 0 ) {
453  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) ) {
454 #ifdef WITHPTHREADS
455  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
456 #endif
457  goto FromNorm;
458  }
459  }
460  ncoef = INCLENG(ncoef);
461  }
462  else if ( d->type == DOLINDEX ) {
463  if ( d->index == 0 ) {
464 #ifdef WITHPTHREADS
465  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
466 #endif
467  goto NormZero;
468  }
469  if ( d->index != NOINDEX ) pind[nind++] = d->index;
470  }
471  else if ( d->type == DOLTERMS ) {
472  t[0] = SUBEXPRESSION;
473  t[4] = AM.dbufnum;
474  if ( t[3] == 0 ) {
475 #ifdef WITHPTHREADS
476  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
477 #endif
478  break;
479  }
480  regval = 2;
481  t = r;
482  while ( t < m ) {
483  if ( *t == DOLLAREXPRESSION ) {
484 #ifdef WITHPTHREADS
485  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
486 #endif
487  d = Dollars + t[2];
488 #ifdef WITHPTHREADS
489  if ( AS.MultiThreaded ) {
490  for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
491  if ( t[2] == ModOptdollars[nummodopt].number ) break;
492  }
493  if ( nummodopt < NumModOptdollars ) {
494  ptype = ModOptdollars[nummodopt].type;
495  if ( ptype == MODLOCAL ) {
496  d = ModOptdollars[nummodopt].dstruct+AT.identity;
497  }
498  else {
499  LOCK(d->pthreadslockread);
500  }
501  }
502  }
503 #endif
504  if ( d->type == DOLTERMS ) {
505  t[0] = SUBEXPRESSION;
506  t[4] = AM.dbufnum;
507  }
508  }
509  t += t[1];
510  }
511 #ifdef WITHPTHREADS
512  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
513 #endif
514  goto RegEnd;
515  }
516  else {
517 #ifdef WITHPTHREADS
518  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
519 #endif
520  MLOCK(ErrorMessageLock);
521  MesPrint("!!!This $ variation has not been implemented yet!!!");
522  MUNLOCK(ErrorMessageLock);
523  goto NormMin;
524  }
525 #ifdef WITHPTHREADS
526  if ( ptype > 0 && ptype != MODLOCAL ) { UNLOCK(d->pthreadslockread); }
527 #endif
528  }
529  else {
530  pnco[nnco++] = t;
531 /*
532  The next statement should be safe as the value is used
533  only by the compiler (ie the master).
534 */
535  AC.lhdollarflag = 1;
536  }
537  break;
538  case DELTA :
539  t += 2;
540  do {
541  if ( *t < 0 ) {
542  if ( *t == SUMMEDIND ) {
543  if ( t[1] < -NMIN4SHIFT ) {
544  k = -t[1]-NMIN4SHIFT;
545  k = ExtraSymbol(k,1,nsym,ppsym,&ncoef);
546  nsym += k;
547  ppsym += (k << 1);
548  }
549  else if ( t[1] == 0 ) goto NormZero;
550  else {
551  if ( t[1] < 0 ) {
552  lnum[0] = -t[1];
553  nnum = -1;
554  }
555  else {
556  lnum[0] = t[1];
557  nnum = 1;
558  }
559  ncoef = REDLENG(ncoef);
560  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) )
561  goto FromNorm;
562  ncoef = INCLENG(ncoef);
563  }
564  t += 2;
565  }
566  else if ( *t == NOINDEX && t[1] == NOINDEX ) t += 2;
567  else if ( *t == EMPTYINDEX && t[1] == EMPTYINDEX ) {
568  *ppdel++ = *t++; *ppdel++ = *t++; ndel += 2;
569  }
570  else
571  if ( t[1] < 0 ) {
572  *ppdot++ = *t++; *ppdot++ = *t++; *ppdot++ = 1; ndot++;
573  }
574  else {
575  *ppvec++ = *t++; *ppvec++ = *t++; nvec += 2;
576  }
577  }
578  else {
579  if ( t[1] < 0 ) {
580  *ppvec++ = t[1]; *ppvec++ = *t; t+=2; nvec += 2;
581  }
582  else { *ppdel++ = *t++; *ppdel++ = *t++; ndel += 2; }
583  }
584  } while ( t < r );
585  break;
586  case FACTORIAL :
587 /*
588  (FACTORIAL,FUNHEAD+2,..,-SNUMBER,number)
589 */
590  if ( t[FUNHEAD] == -SNUMBER && t[1] == FUNHEAD+2
591  && t[FUNHEAD+1] >= 0 ) {
592  if ( Factorial(BHEAD t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
593  goto FromNorm;
594 MulIn:
595  ncoef = REDLENG(ncoef);
596  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) ) goto FromNorm;
597  ncoef = INCLENG(ncoef);
598  }
599  else pcom[ncom++] = t;
600  break;
601  case BERNOULLIFUNCTION :
602 /*
603  (BERNOULLIFUNCTION,FUNHEAD+2,..,-SNUMBER,number)
604 */
605  if ( ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] >= 0 )
606  && ( t[1] == FUNHEAD+2 || ( t[1] == FUNHEAD+4 &&
607  t[FUNHEAD+2] == -SNUMBER && ABS(t[FUNHEAD+3]) == 1 ) ) ) {
608  WORD inum, mnum;
609  if ( Bernoulli(t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
610  goto FromNorm;
611  if ( nnum == 0 ) goto NormZero;
612  inum = nnum; if ( inum < 0 ) inum = -inum;
613  inum--; inum /= 2;
614  mnum = inum;
615  while ( lnum[mnum-1] == 0 ) mnum--;
616  if ( nnum < 0 ) mnum = -mnum;
617  ncoef = REDLENG(ncoef);
618  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,mnum) ) goto FromNorm;
619  mnum = inum;
620  while ( lnum[inum+mnum-1] == 0 ) mnum--;
621  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)(lnum+inum),mnum) ) goto FromNorm;
622  ncoef = INCLENG(ncoef);
623  if ( t[1] == FUNHEAD+4 && t[FUNHEAD+1] == 1
624  && t[FUNHEAD+3] == -1 ) ncoef = -ncoef;
625  }
626  else pcom[ncom++] = t;
627  break;
628  case NUMARGSFUN:
629 /*
630  Numerical function giving the number of arguments.
631 */
632  k = 0;
633  t += FUNHEAD;
634  while ( t < r ) {
635  k++;
636  NEXTARG(t);
637  }
638  if ( k == 0 ) goto NormZero;
639  *((UWORD *)lnum) = k;
640  nnum = 1;
641  goto MulIn;
642  case NUMFACTORS:
643 /*
644  Numerical function giving the number of factors in an expression.
645 */
646  t += FUNHEAD;
647  if ( *t == -EXPRESSION ) {
648  k = AS.OldNumFactors[t[1]];
649  }
650  else if ( *t == -DOLLAREXPRESSION ) {
651  k = Dollars[t[1]].nfactors;
652  }
653  else {
654  pcom[ncom++] = t;
655  break;
656  }
657  if ( k == 0 ) goto NormZero;
658  *((UWORD *)lnum) = k;
659  nnum = 1;
660  goto MulIn;
661  case NUMTERMSFUN:
662 /*
663  Numerical function giving the number of terms in the single argument.
664 */
665  if ( t[FUNHEAD] < 0 ) {
666  if ( t[FUNHEAD] <= -FUNCTION && t[1] == FUNHEAD+1 ) break;
667  if ( t[FUNHEAD] > -FUNCTION && t[1] == FUNHEAD+2 ) {
668  if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] == 0 ) goto NormZero;
669  break;
670  }
671  pcom[ncom++] = t;
672  break;
673  }
674  if ( t[FUNHEAD] > 0 && t[FUNHEAD] == t[1]-FUNHEAD ) {
675  k = 0;
676  t += FUNHEAD+ARGHEAD;
677  while ( t < r ) {
678  k++;
679  t += *t;
680  }
681  if ( k == 0 ) goto NormZero;
682  *((UWORD *)lnum) = k;
683  nnum = 1;
684  goto MulIn;
685  }
686  else pcom[ncom++] = t;
687  break;
688  case COUNTFUNCTION:
689  if ( AN.cTerm ) {
690  k = CountFun(AN.cTerm,t);
691  if ( k == 0 ) goto NormZero;
692  if ( k > 0 ) { *((UWORD *)lnum) = k; nnum = 1; }
693  else { *((UWORD *)lnum) = -k; nnum = -1; }
694  goto MulIn;
695  }
696  break;
697  case MAKERATIONAL:
698  if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+2] == -SNUMBER
699  && t[1] == FUNHEAD+4 && t[FUNHEAD+3] > 1 ) {
700  WORD x1[2], sgn;
701  if ( t[FUNHEAD+1] == 0 ) goto NormZero;
702  if ( t[FUNHEAD+1] < 0 ) { t[FUNHEAD+1] = -t[FUNHEAD+1]; sgn = -1; }
703  else sgn = 1;
704  if ( MakeRational(t[FUNHEAD+1],t[FUNHEAD+3],x1,x1+1) ) {
705  static int warnflag = 1;
706  if ( warnflag ) {
707  MesPrint("%w Warning: fraction could not be reconstructed in MakeRational_");
708  warnflag = 0;
709  }
710  x1[0] = t[FUNHEAD+1]; x1[1] = 1;
711  }
712  if ( sgn < 0 ) { t[FUNHEAD+1] = -t[FUNHEAD+1]; x1[0] = -x1[0]; }
713  if ( x1[0] < 0 ) { sgn = -1; x1[0] = -x1[0]; }
714  else sgn = 1;
715  ncoef = REDLENG(ncoef);
716  if ( MulRat(BHEAD (UWORD *)AT.n_coef,ncoef,(UWORD *)x1,sgn,
717  (UWORD *)AT.n_coef,&ncoef) ) goto FromNorm;
718  ncoef = INCLENG(ncoef);
719  }
720  else {
721  WORD narg = 0, *tt, *ttstop, *arg1 = 0, *arg2 = 0;
722  UWORD *x1, *x2, *xx;
723  WORD nx1,nx2,nxx;
724  ttstop = t + t[1]; tt = t+FUNHEAD;
725  while ( tt < ttstop ) {
726  narg++;
727  if ( narg == 1 ) arg1 = tt;
728  else arg2 = tt;
729  NEXTARG(tt);
730  }
731  if ( narg != 2 ) goto defaultcase;
732  if ( *arg2 == -SNUMBER && arg2[1] <= 1 ) goto defaultcase;
733  else if ( *arg2 > 0 && ttstop[-1] < 0 ) goto defaultcase;
734  x1 = NumberMalloc("Norm-MakeRational");
735  if ( *arg1 == -SNUMBER ) {
736  if ( arg1[1] == 0 ) {
737  NumberFree(x1,"Norm-MakeRational");
738  goto NormZero;
739  }
740  if ( arg1[1] < 0 ) { x1[0] = -arg1[1]; nx1 = -1; }
741  else { x1[0] = arg1[1]; nx1 = 1; }
742  }
743  else if ( *arg1 > 0 ) {
744  WORD *tc;
745  nx1 = (ABS(arg2[-1])-1)/2;
746  tc = arg1+ARGHEAD+1+nx1;
747  if ( tc[0] != 1 ) {
748  NumberFree(x1,"Norm-MakeRational");
749  goto defaultcase;
750  }
751  for ( i = 1; i < nx1; i++ ) if ( tc[i] != 0 ) {
752  NumberFree(x1,"Norm-MakeRational");
753  goto defaultcase;
754  }
755  tc = arg1+ARGHEAD+1;
756  for ( i = 0; i < nx1; i++ ) x1[i] = tc[i];
757  if ( arg2[-1] < 0 ) nx1 = -nx1;
758  }
759  else {
760  NumberFree(x1,"Norm-MakeRational");
761  goto defaultcase;
762  }
763  x2 = NumberMalloc("Norm-MakeRational");
764  if ( *arg2 == -SNUMBER ) {
765  if ( arg2[1] <= 1 ) {
766  NumberFree(x2,"Norm-MakeRational");
767  NumberFree(x1,"Norm-MakeRational");
768  goto defaultcase;
769  }
770  else { x2[0] = arg2[1]; nx2 = 1; }
771  }
772  else if ( *arg2 > 0 ) {
773  WORD *tc;
774  nx2 = (ttstop[-1]-1)/2;
775  tc = arg2+ARGHEAD+1+nx2;
776  if ( tc[0] != 1 ) {
777  NumberFree(x2,"Norm-MakeRational");
778  NumberFree(x1,"Norm-MakeRational");
779  goto defaultcase;
780  }
781  for ( i = 1; i < nx2; i++ ) if ( tc[i] != 0 ) {
782  NumberFree(x2,"Norm-MakeRational");
783  NumberFree(x1,"Norm-MakeRational");
784  goto defaultcase;
785  }
786  tc = arg2+ARGHEAD+1;
787  for ( i = 0; i < nx2; i++ ) x2[i] = tc[i];
788  }
789  else {
790  NumberFree(x2,"Norm-MakeRational");
791  NumberFree(x1,"Norm-MakeRational");
792  goto defaultcase;
793  }
794  if ( BigLong(x1,ABS(nx1),x2,nx2) >= 0 ) {
795  UWORD *x3 = NumberMalloc("Norm-MakeRational");
796  UWORD *x4 = NumberMalloc("Norm-MakeRational");
797  WORD nx3, nx4;
798  DivLong(x1,nx1,x2,nx2,x3,&nx3,x4,&nx4);
799  for ( i = 0; i < ABS(nx4); i++ ) x1[i] = x4[i];
800  nx1 = nx4;
801  NumberFree(x4,"Norm-MakeRational");
802  NumberFree(x3,"Norm-MakeRational");
803  }
804  xx = (UWORD *)(TermMalloc("Norm-MakeRational"));
805  if ( MakeLongRational(BHEAD x1,nx1,x2,nx2,xx,&nxx) ) {
806  static int warnflag = 1;
807  if ( warnflag ) {
808  MesPrint("%w Warning: fraction could not be reconstructed in MakeRational_");
809  warnflag = 0;
810  }
811  ncoef = REDLENG(ncoef);
812  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,x1,nx1) )
813  goto FromNorm;
814  }
815  else {
816  ncoef = REDLENG(ncoef);
817  if ( MulRat(BHEAD (UWORD *)AT.n_coef,ncoef,xx,nxx,
818  (UWORD *)AT.n_coef,&ncoef) ) goto FromNorm;
819  }
820  ncoef = INCLENG(ncoef);
821  TermFree(xx,"Norm-MakeRational");
822  NumberFree(x2,"Norm-MakeRational");
823  NumberFree(x1,"Norm-MakeRational");
824  }
825  break;
826  case TERMFUNCTION:
827  if ( t[1] == FUNHEAD && AN.cTerm ) {
828  ANsr = r; ANsm = m; ANsc = AN.cTerm;
829  AN.cTerm = 0;
830  t = ANsc + 1;
831  m = ANsc + *ANsc;
832  ncoef = REDLENG(ncoef);
833  nnum = REDLENG(m[-1]);
834  m -= ABS(m[-1]);
835  if ( MulRat(BHEAD (UWORD *)AT.n_coef,ncoef,(UWORD *)m,nnum,
836  (UWORD *)AT.n_coef,&ncoef) ) goto FromNorm;
837  ncoef = INCLENG(ncoef);
838  r = t;
839  }
840  break;
841  case FIRSTBRACKET:
842  if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
843  if ( GetFirstBracket(termout,t[FUNHEAD+1]) < 0 ) goto FromNorm;
844  if ( *termout == 0 ) goto NormZero;
845  if ( *termout > 4 ) {
846  WORD *r1, *r2, *r3;
847  while ( r < m ) *t++ = *r++;
848  r1 = term + *term;
849  r2 = termout + *termout; r2 -= ABS(r2[-1]);
850  while ( r < r1 ) *r2++ = *r++;
851  r3 = termout + 1;
852  while ( r3 < r2 ) *t++ = *r3++; *term = t - term;
853  if ( AT.WorkPointer > term && AT.WorkPointer < t )
854  AT.WorkPointer = t;
855  goto Restart;
856  }
857  }
858  break;
859  case FIRSTTERM:
860  case CONTENTTERM:
861  if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
862  {
863  EXPRESSIONS e = Expressions+t[FUNHEAD+1];
864  POSITION oldondisk = AS.OldOnFile[t[FUNHEAD+1]];
865  if ( e->replace == NEWLYDEFINEDEXPRESSION ) {
866  AS.OldOnFile[t[FUNHEAD+1]] = e->onfile;
867  }
868  if ( *t == FIRSTTERM ) {
869  if ( GetFirstTerm(termout,t[FUNHEAD+1]) < 0 ) goto FromNorm;
870  }
871  else if ( *t == CONTENTTERM ) {
872  if ( GetContent(termout,t[FUNHEAD+1]) < 0 ) goto FromNorm;
873  }
874  AS.OldOnFile[t[FUNHEAD+1]] = oldondisk;
875  if ( *termout == 0 ) goto NormZero;
876  }
877 PasteIn:;
878  {
879  WORD *r1, *r2, *r3, *r4, *r5, nr1, *rterm;
880  r2 = termout + *termout; lnum = r2 - ABS(r2[-1]);
881  nnum = REDLENG(r2[-1]);
882 
883  r1 = term + *term; r3 = r1 - ABS(r1[-1]);
884  nr1 = REDLENG(r1[-1]);
885  if ( Mully(BHEAD (UWORD *)lnum,&nnum,(UWORD *)r3,nr1) ) goto FromNorm;
886  nnum = INCLENG(nnum); nr1 = ABS(nnum); lnum[nr1-1] = nnum;
887  rterm = TermMalloc("FirstTerm/ContentTerm");
888  r4 = rterm+1; r5 = term+1; while ( r5 < t ) *r4++ = *r5++;
889  r5 = termout+1; while ( r5 < lnum ) *r4++ = *r5++;
890  r5 = r; while ( r5 < r3 ) *r4++ = *r5++;
891  r5 = lnum; NCOPY(r4,r5,nr1);
892  *rterm = r4-rterm;
893  nr1 = *rterm; r1 = term; r2 = rterm; NCOPY(r1,r2,nr1);
894  TermFree(rterm,"FirstTerm/ContentTerm");
895  if ( AT.WorkPointer > term && AT.WorkPointer < r1 )
896  AT.WorkPointer = r1;
897  goto Restart;
898  }
899  }
900  else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
901  DOLLARS d = Dollars + t[FUNHEAD+1], newd = 0;
902  int idol, ido;
903 #ifdef WITHPTHREADS
904  int nummodopt, dtype = -1;
905  if ( AS.MultiThreaded && ( AC.mparallelflag == PARALLELFLAG ) ) {
906  for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
907  if ( t[FUNHEAD+1] == ModOptdollars[nummodopt].number ) break;
908  }
909  if ( nummodopt < NumModOptdollars ) {
910  dtype = ModOptdollars[nummodopt].type;
911  if ( dtype == MODLOCAL ) {
912  d = ModOptdollars[nummodopt].dstruct+AT.identity;
913  }
914  }
915  }
916 #endif
917  if ( d->where && ( d->type == DOLTERMS || d->type == DOLNUMBER ) ) {
918  newd = d;
919  }
920  else {
921  if ( ( newd = DolToTerms(BHEAD t[FUNHEAD+1]) ) == 0 )
922  goto NormZero;
923  }
924  if ( newd->where[0] == 0 ) {
925  M_free(newd,"Copy of dollar variable");
926  goto NormZero;
927  }
928  if ( *t == FIRSTTERM ) {
929  idol = newd->where[0];
930  for ( ido = 0; ido < idol; ido++ ) termout[ido] = newd->where[ido];
931  }
932  else if ( *t == CONTENTTERM ) {
933  WORD *tterm;
934  tterm = newd->where;
935  idol = tterm[0];
936  for ( ido = 0; ido < idol; ido++ ) termout[ido] = tterm[ido];
937  tterm += *tterm;
938  while ( *tterm ) {
939  if ( ContentMerge(BHEAD termout,tterm) < 0 ) goto FromNorm;
940  tterm += *tterm;
941  }
942  }
943  if ( newd != d ) {
944  if ( newd->factors ) M_free(newd->factors,"Dollar factors");
945  M_free(newd,"Copy of dollar variable");
946  newd = 0;
947  }
948  goto PasteIn;
949  }
950  break;
951  case TERMSINEXPR:
952  { LONG x;
953  if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
954  x = TermsInExpression(t[FUNHEAD+1]);
955 multermnum: if ( x == 0 ) goto NormZero;
956  if ( x < 0 ) {
957  x = -x;
958  if ( x > (LONG)WORDMASK ) { lnum[0] = x & WORDMASK;
959  lnum[1] = x >> BITSINWORD; nnum = -2;
960  }
961  else { lnum[0] = x; nnum = -1; }
962  }
963  else if ( x > (LONG)WORDMASK ) {
964  lnum[0] = x & WORDMASK;
965  lnum[1] = x >> BITSINWORD;
966  nnum = 2;
967  }
968  else { lnum[0] = x; nnum = 1; }
969  ncoef = REDLENG(ncoef);
970  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) )
971  goto FromNorm;
972  ncoef = INCLENG(ncoef);
973  }
974  else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
975  x = TermsInDollar(t[FUNHEAD+1]);
976  goto multermnum;
977  }
978  else { pcom[ncom++] = t; }
979  }
980  break;
981  case MATCHFUNCTION:
982  case PATTERNFUNCTION:
983  break;
984  case BINOMIAL:
985 /*
986  Binomial function for internal use for the moment.
987  The routine in reken.c should be more efficient.
988 */
989  if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
990  && t[FUNHEAD+1] >= 0 && t[FUNHEAD+2] == -SNUMBER
991  && t[FUNHEAD+3] >= 0 && t[FUNHEAD+1] >= t[FUNHEAD+3] ) {
992  if ( t[FUNHEAD+1] > t[FUNHEAD+3] ) {
993  if ( GetBinom((UWORD *)lnum,&nnum,
994  t[FUNHEAD+1],t[FUNHEAD+3]) ) goto FromNorm;
995  if ( nnum == 0 ) goto NormZero;
996  goto MulIn;
997  }
998  }
999  else pcom[ncom++] = t;
1000  break;
1001  case SIGNFUN:
1002 /*
1003  Numerical function giving (-1)^arg
1004 */
1005  if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1006  if ( ( t[FUNHEAD+1] & 1 ) != 0 ) ncoef = -ncoef;
1007  }
1008  else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] )
1009  && ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) ) {
1010  UWORD *numer1,*denom1;
1011  WORD nsize = abs(t[t[1]-1]), nnsize, isize;
1012  nnsize = (nsize-1)/2;
1013  numer1 = (UWORD *)(t + FUNHEAD+ARGHEAD+1);
1014  denom1 = numer1 + nnsize;
1015  for ( isize = 1; isize < nnsize; isize++ ) {
1016  if ( denom1[isize] ) break;
1017  }
1018  if ( ( denom1[0] != 1 ) || isize < nnsize ) {
1019  pcom[ncom++] = t;
1020  }
1021  else {
1022  if ( ( numer1[0] & 1 ) != 0 ) ncoef = -ncoef;
1023  }
1024  }
1025  else {
1026  goto doflags;
1027 /* pcom[ncom++] = t; */
1028  }
1029  break;
1030  case SIGFUNCTION:
1031 /*
1032  Numerical function giving the sign of the numerical argument
1033  The sign of zero is 1.
1034  If there are roots of unity they are part of the sign.
1035 */
1036  if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1037  if ( t[FUNHEAD+1] < 0 ) ncoef = -ncoef;
1038  }
1039  else if ( ( t[1] == FUNHEAD+2 ) && ( t[FUNHEAD] == -SYMBOL )
1040  && ( ( t[FUNHEAD+1] <= NumSymbols && t[FUNHEAD+1] > -MAXPOWER )
1041  && ( symbols[t[FUNHEAD+1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) ) {
1042  k = t[FUNHEAD+1];
1043  from = m;
1044  i = nsym;
1045  m = ppsym;
1046  if ( i > 0 ) do {
1047  m -= 2;
1048  if ( k == *m ) {
1049  m++;
1050  *m = *m + 1;
1051  *m %= symbols[k].maxpower;
1052  if ( ( symbols[k].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
1053  if ( ( ( symbols[k].maxpower & 1 ) == 0 ) &&
1054  ( *m >= symbols[k].maxpower/2 ) ) {
1055  *m -= symbols[k].maxpower/2; ncoef = -ncoef;
1056  }
1057  }
1058  if ( !*m ) {
1059  m--;
1060  while ( i < nsym )
1061  { *m = m[2]; m++; *m = m[2]; m++; i++; }
1062  ppsym -= 2;
1063  nsym--;
1064  }
1065  goto sigDoneSymbol;
1066  }
1067  } while ( k < *m && --i > 0 );
1068  m = ppsym;
1069  while ( i < nsym )
1070  { m--; m[2] = *m; m--; m[2] = *m; i++; }
1071  *m++ = k;
1072  *m = 1;
1073  ppsym += 2;
1074  nsym++;
1075 sigDoneSymbol:;
1076  m = from;
1077  }
1078  else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] ) ) {
1079  if ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) {
1080  if ( t[t[1]-1] < 0 ) ncoef = -ncoef;
1081  }
1082 /*
1083  Now we should fish out the roots of unity
1084 */
1085  else if ( ( t[FUNHEAD+ARGHEAD]+FUNHEAD+ARGHEAD == t[1] )
1086  && ( t[FUNHEAD+ARGHEAD+1] == SYMBOL ) ) {
1087  WORD *ts = t + FUNHEAD+ARGHEAD+3;
1088  WORD its = ts[-1]-2;
1089  while ( its > 0 ) {
1090  if ( ( *ts != 0 ) && ( ( *ts > NumSymbols || *ts <= -MAXPOWER )
1091  || ( symbols[*ts].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY ) ) {
1092  goto signogood;
1093  }
1094  ts += 2; its -= 2;
1095  }
1096 /*
1097  Now we have only roots of unity which should be
1098  registered in the list of sysmbols.
1099 */
1100  if ( t[t[1]-1] < 0 ) ncoef = -ncoef;
1101  ts = t + FUNHEAD+ARGHEAD+3;
1102  its = ts[-1]-2;
1103  from = m;
1104  while ( its > 0 ) {
1105  i = nsym;
1106  m = ppsym;
1107  if ( i > 0 ) do {
1108  m -= 2;
1109  if ( *ts == *m ) {
1110  ts++; m++;
1111  *m += *ts;
1112  if ( ( ts[-1] <= NumSymbols && ts[-1] > -MAXPOWER ) &&
1113  ( symbols[ts[-1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
1114  *m %= symbols[ts[-1]].maxpower;
1115  if ( *m < 0 ) *m += symbols[ts[-1]].maxpower;
1116  if ( ( symbols[ts[-1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
1117  if ( ( ( symbols[ts[-1]].maxpower & 1 ) == 0 ) &&
1118  ( *m >= symbols[ts[-1]].maxpower/2 ) ) {
1119  *m -= symbols[ts[-1]].maxpower/2; ncoef = -ncoef;
1120  }
1121  }
1122  }
1123  if ( !*m ) {
1124  m--;
1125  while ( i < nsym )
1126  { *m = m[2]; m++; *m = m[2]; m++; i++; }
1127  ppsym -= 2;
1128  nsym--;
1129  }
1130  ts++; its -= 2;
1131  goto sigNextSymbol;
1132  }
1133  } while ( *ts < *m && --i > 0 );
1134  m = ppsym;
1135  while ( i < nsym )
1136  { m--; m[2] = *m; m--; m[2] = *m; i++; }
1137  *m++ = *ts++;
1138  *m = *ts++;
1139  ppsym += 2;
1140  nsym++; its -= 2;
1141 sigNextSymbol:;
1142  }
1143  m = from;
1144  }
1145  else {
1146 signogood: pcom[ncom++] = t;
1147  }
1148  }
1149  else pcom[ncom++] = t;
1150  break;
1151  case ABSFUNCTION:
1152 /*
1153  Numerical function giving the absolute value of the
1154  numerical argument. Or roots of unity.
1155 */
1156  if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1157  k = t[FUNHEAD+1];
1158  if ( k < 0 ) k = -k;
1159  if ( k == 0 ) goto NormZero;
1160  *((UWORD *)lnum) = k; nnum = 1;
1161  goto MulIn;
1162 
1163  }
1164  else if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SYMBOL ) {
1165  k = t[FUNHEAD+1];
1166  if ( ( k > NumSymbols || k <= -MAXPOWER )
1167  || ( symbols[k].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
1168  goto absnogood;
1169  }
1170  else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] )
1171  && ( t[1] == FUNHEAD+ARGHEAD+t[FUNHEAD+ARGHEAD] ) ) {
1172  if ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) {
1173  WORD *ts;
1174 absnosymbols: ts = t + t[1] -1;
1175  ncoef = REDLENG(ncoef);
1176  nnum = REDLENG(*ts);
1177  if ( nnum < 0 ) nnum = -nnum;
1178  if ( MulRat(BHEAD (UWORD *)AT.n_coef,ncoef,
1179  (UWORD *)(ts-ABS(*ts)+1),nnum,
1180  (UWORD *)AT.n_coef,&ncoef) ) goto FromNorm;
1181  ncoef = INCLENG(ncoef);
1182  }
1183 /*
1184  Now get rid of the roots of unity. This includes i_
1185 */
1186  else if ( t[FUNHEAD+ARGHEAD+1] == SYMBOL ) {
1187  WORD *ts = t+FUNHEAD+ARGHEAD+1;
1188  WORD its = ts[1] - 2;
1189  ts += 2;
1190  while ( its > 0 ) {
1191  if ( *ts == 0 ) { }
1192  else if ( ( *ts > NumSymbols || *ts <= -MAXPOWER )
1193  || ( symbols[*ts].complex & VARTYPEROOTOFUNITY )
1194  != VARTYPEROOTOFUNITY ) goto absnogood;
1195  its -= 2; ts += 2;
1196  }
1197  goto absnosymbols;
1198  }
1199  else {
1200 absnogood: pcom[ncom++] = t;
1201  }
1202  }
1203  else pcom[ncom++] = t;
1204  break;
1205  case MODFUNCTION:
1206  case MOD2FUNCTION:
1207 /*
1208  Mod function. Does work if two arguments and the
1209  second argument is a positive short number
1210 */
1211  if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1212  && t[FUNHEAD+2] == -SNUMBER && t[FUNHEAD+3] > 1 ) {
1213  WORD tmod;
1214  tmod = (t[FUNHEAD+1]%t[FUNHEAD+3]);
1215  if ( tmod < 0 ) tmod += t[FUNHEAD+3];
1216  if ( *t == MOD2FUNCTION && tmod > t[FUNHEAD+3]/2 )
1217  tmod -= t[FUNHEAD+3];
1218  if ( tmod < 0 ) {
1219  *((UWORD *)lnum) = -tmod;
1220  nnum = -1;
1221  }
1222  else if ( tmod > 0 ) {
1223  *((UWORD *)lnum) = tmod;
1224  nnum = 1;
1225  }
1226  else goto NormZero;
1227  goto MulIn;
1228  }
1229  else if ( t[1] > t[FUNHEAD+2] && t[FUNHEAD] > 0
1230  && t[FUNHEAD+t[FUNHEAD]] == -SNUMBER
1231  && t[FUNHEAD+t[FUNHEAD]+1] > 1
1232  && t[1] == FUNHEAD+2+t[FUNHEAD] ) {
1233  WORD *ttt = t+FUNHEAD, iii;
1234  iii = ttt[*ttt-1];
1235  if ( *ttt == ttt[ARGHEAD]+ARGHEAD &&
1236  ttt[ARGHEAD] == ABS(iii)+1 ) {
1237  WORD ncmod = 1;
1238  WORD cmod = ttt[*ttt+1];
1239  iii = REDLENG(iii);
1240  if ( *t == MODFUNCTION ) {
1241  if ( TakeModulus((UWORD *)(ttt+ARGHEAD+1)
1242  ,&iii,(UWORD *)(&cmod),ncmod,UNPACK|NOINVERSES) )
1243  goto FromNorm;
1244  }
1245  else {
1246  if ( TakeModulus((UWORD *)(ttt+ARGHEAD+1)
1247  ,&iii,(UWORD *)(&cmod),ncmod,UNPACK|POSNEG|NOINVERSES) )
1248  goto FromNorm;
1249  }
1250  if ( *t == MOD2FUNCTION && ttt[ARGHEAD+1] > cmod/2 && iii > 0 ) {
1251  ttt[ARGHEAD+1] -= cmod;
1252  }
1253  if ( ttt[ARGHEAD+1] < 0 ) {
1254  *((UWORD *)lnum) = -ttt[ARGHEAD+1];
1255  nnum = -1;
1256  }
1257  else if ( ttt[ARGHEAD+1] > 0 ) {
1258  *((UWORD *)lnum) = ttt[ARGHEAD+1];
1259  nnum = 1;
1260  }
1261  else goto NormZero;
1262  goto MulIn;
1263  }
1264  }
1265  else if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1266  *((UWORD *)lnum) = t[FUNHEAD+1];
1267  if ( *lnum == 0 ) goto NormZero;
1268  nnum = 1;
1269  goto MulIn;
1270  }
1271  else if ( ( ( t[FUNHEAD] < 0 ) && ( t[FUNHEAD] == -SNUMBER )
1272  && ( t[1] >= ( FUNHEAD+6+ARGHEAD ) )
1273  && ( t[FUNHEAD+2] >= 4+ARGHEAD )
1274  && ( t[t[1]-1] == t[FUNHEAD+2+ARGHEAD]-1 ) ) ||
1275  ( ( t[FUNHEAD] > 0 )
1276  && ( t[FUNHEAD]-ARGHEAD-1 == ABS(t[FUNHEAD+t[FUNHEAD]-1]) )
1277  && ( t[FUNHEAD+t[FUNHEAD]]-ARGHEAD-1 == t[t[1]-1] ) ) ) {
1278 /*
1279  Check that the last (long) number is integer
1280 */
1281  WORD *ttt = t + t[1], iii, iii1;
1282  UWORD coefbuf[2], *coef2, ncoef2;
1283  iii = (ttt[-1]-1)/2;
1284  ttt -= iii;
1285  if ( ttt[-1] != 1 ) {
1286 exitfromhere:
1287  pcom[ncom++] = t;
1288  break;
1289  }
1290  iii--;
1291  for ( iii1 = 0; iii1 < iii; iii1++ ) {
1292  if ( ttt[iii1] != 0 ) goto exitfromhere;
1293  }
1294 /*
1295  Now we have a hit!
1296  The first argument will be put in lnum.
1297  It will be a rational.
1298  The second argument will be a long integer in coef2.
1299 */
1300  ttt = t + FUNHEAD;
1301  if ( *ttt < 0 ) {
1302  if ( ttt[1] < 0 ) {
1303  nnum = -1; lnum[0] = -ttt[1]; lnum[1] = 1;
1304  }
1305  else {
1306  nnum = 1; lnum[0] = ttt[1]; lnum[1] = 1;
1307  }
1308  }
1309  else {
1310  nnum = ABS(ttt[ttt[0]-1] - 1);
1311  for ( iii = 0; iii < nnum; iii++ ) {
1312  lnum[iii] = ttt[ARGHEAD+1+iii];
1313  }
1314  nnum = nnum/2;
1315  if ( ttt[ttt[0]-1] < 0 ) nnum = -nnum;
1316  }
1317  NEXTARG(ttt);
1318  if ( *ttt < 0 ) {
1319  coef2 = coefbuf;
1320  ncoef2 = 3; *coef2 = (UWORD)(ttt[1]);
1321  coef2[1] = 1;
1322  }
1323  else {
1324  coef2 = (UWORD *)(ttt+ARGHEAD+1);
1325  ncoef2 = (ttt[ttt[0]-1]-1)/2;
1326  }
1327  if ( TakeModulus((UWORD *)lnum,&nnum,(UWORD *)coef2,ncoef2,
1328  UNPACK|NOINVERSES|FROMFUNCTION) ) {
1329  goto FromNorm;
1330  }
1331  if ( *t == MOD2FUNCTION && nnum > 0 ) {
1332  UWORD *coef3 = NumberMalloc("Mod2Function"), two = 2;
1333  WORD ncoef3;
1334  if ( MulLong((UWORD *)lnum,nnum,&two,1,coef3,&ncoef3) )
1335  goto FromNorm;
1336  if ( BigLong(coef3,ncoef3,(UWORD *)coef2,ncoef2) > 0 ) {
1337  nnum = -nnum;
1338  AddLong((UWORD *)lnum,nnum,(UWORD *)coef2,ncoef2
1339  ,(UWORD *)lnum,&nnum);
1340  nnum = -nnum;
1341  }
1342  NumberFree(coef3,"Mod2Function");
1343  }
1344 /*
1345  Do we have to pack? No, because the answer is not a fraction
1346 */
1347  ncoef = REDLENG(ncoef);
1348  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) )
1349  goto FromNorm;
1350  ncoef = INCLENG(ncoef);
1351  }
1352  else pcom[ncom++] = t;
1353  break;
1354  case EXTEUCLIDEAN:
1355  {
1356  WORD argcount = 0, *tc, *ts, xc, xs, *tcc;
1357  UWORD *Num1, *Num2, *Num3, *Num4;
1358  WORD size1, size2, size3, size4, space;
1359  tc = t+FUNHEAD; ts = t + t[1];
1360  while ( argcount < 3 && tc < ts ) { NEXTARG(tc); argcount++; }
1361  if ( argcount != 2 ) goto defaultcase;
1362  if ( t[FUNHEAD] == -SNUMBER ) {
1363  if ( t[FUNHEAD+1] <= 1 ) goto defaultcase;
1364  if ( t[FUNHEAD+2] == -SNUMBER ) {
1365  if ( t[FUNHEAD+3] <= 1 ) goto defaultcase;
1366  Num2 = NumberMalloc("modinverses");
1367  *Num2 = t[FUNHEAD+3]; size2 = 1;
1368  }
1369  else {
1370  if ( ts[-1] < 0 ) goto defaultcase;
1371  if ( ts[-1] != t[FUNHEAD+2]-ARGHEAD-1 ) goto defaultcase;
1372  xs = (ts[-1]-1)/2;
1373  tcc = ts-xs-1;
1374  if ( *tcc != 1 ) goto defaultcase;
1375  for ( i = 1; i < xs; i++ ) {
1376  if ( tcc[i] != 0 ) goto defaultcase;
1377  }
1378  Num2 = NumberMalloc("modinverses");
1379  size2 = xs;
1380  for ( i = 0; i < xs; i++ ) Num2[i] = t[FUNHEAD+ARGHEAD+3+i];
1381  }
1382  Num1 = NumberMalloc("modinverses");
1383  *Num1 = t[FUNHEAD+1]; size1 = 1;
1384  }
1385  else {
1386  tc = t + FUNHEAD + t[FUNHEAD];
1387  if ( tc[-1] < 0 ) goto defaultcase;
1388  if ( tc[-1] != t[FUNHEAD]-ARGHEAD-1 ) goto defaultcase;
1389  xc = (tc[-1]-1)/2;
1390  tcc = tc-xc-1;
1391  if ( *tcc != 1 ) goto defaultcase;
1392  for ( i = 1; i < xc; i++ ) {
1393  if ( tcc[i] != 0 ) goto defaultcase;
1394  }
1395  if ( *tc == -SNUMBER ) {
1396  if ( tc[1] <= 1 ) goto defaultcase;
1397  Num2 = NumberMalloc("modinverses");
1398  *Num2 = tc[1]; size2 = 1;
1399  }
1400  else {
1401  if ( ts[-1] < 0 ) goto defaultcase;
1402  if ( ts[-1] != t[FUNHEAD+2]-ARGHEAD-1 ) goto defaultcase;
1403  xs = (ts[-1]-1)/2;
1404  tcc = ts-xs-1;
1405  if ( *tcc != 1 ) goto defaultcase;
1406  for ( i = 1; i < xs; i++ ) {
1407  if ( tcc[i] != 0 ) goto defaultcase;
1408  }
1409  Num2 = NumberMalloc("modinverses");
1410  size2 = xs;
1411  for ( i = 0; i < xs; i++ ) Num2[i] = tc[ARGHEAD+1+i];
1412  }
1413  Num1 = NumberMalloc("modinverses");
1414  size1 = xc;
1415  for ( i = 0; i < xc; i++ ) Num1[i] = t[FUNHEAD+ARGHEAD+1+i];
1416  }
1417  Num3 = NumberMalloc("modinverses");
1418  Num4 = NumberMalloc("modinverses");
1419  GetLongModInverses(BHEAD Num1,size1,Num2,size2
1420  ,Num3,&size3,Num4,&size4);
1421 /*
1422  Now we have to compose the answer. This needs more space
1423  and hence we have to put this inside the term.
1424  Compute first how much extra space we need.
1425  Then move the trailing part of the term upwards.
1426  Do not forget relevant pointers!!! (r, m, termout, AT.WorkPointer)
1427 */
1428  space = 0;
1429  if ( ( size3 == 1 || size3 == -1 ) && (*Num3&TOPBITONLY) == 0 ) space += 2;
1430  else space += ARGHEAD + 2*ABS(size3) + 2;
1431  if ( ( size4 == 1 || size4 == -1 ) && (*Num4&TOPBITONLY) == 0 ) space += 2;
1432  else space += ARGHEAD + 2*ABS(size4) + 2;
1433  tt = term + *term; u = tt + space;
1434  while ( tt >= ts ) *--u = *--tt;
1435  m += space; r += space;
1436  *term += space;
1437  t[1] += space;
1438  if ( ( size3 == 1 || size3 == -1 ) && (*Num3&TOPBITONLY) == 0 ) {
1439  *ts++ = -SNUMBER; *ts = (WORD)(*Num3);
1440  if ( size3 < 0 ) *ts = -*ts;
1441  ts++;
1442  }
1443  else {
1444  *ts++ = 2*ABS(size3)+ARGHEAD+2;
1445  *ts++ = 0; FILLARG(ts)
1446  *ts++ = 2*ABS(size3)+1;
1447  for ( i = 0; i < ABS(size3); i++ ) *ts++ = Num3[i];
1448  *ts++ = 1;
1449  for ( i = 1; i < ABS(size3); i++ ) *ts++ = 0;
1450  if ( size3 < 0 ) *ts++ = 2*size3-1;
1451  else *ts++ = 2*size3+1;
1452  }
1453  if ( ( size4 == 1 || size4 == -1 ) && (*Num4&TOPBITONLY) == 0 ) {
1454  *ts++ = -SNUMBER; *ts = *Num4;
1455  if ( size4 < 0 ) *ts = -*ts;
1456  ts++;
1457  }
1458  else {
1459  *ts++ = 2*ABS(size4)+ARGHEAD+2;
1460  *ts++ = 0; FILLARG(ts)
1461  *ts++ = 2*ABS(size4)+2;
1462  for ( i = 0; i < ABS(size4); i++ ) *ts++ = Num4[i];
1463  *ts++ = 1;
1464  for ( i = 1; i < ABS(size4); i++ ) *ts++ = 0;
1465  if ( size4 < 0 ) *ts++ = 2*size4-1;
1466  else *ts++ = 2*size4+1;
1467  }
1468  NumberFree(Num4,"modinverses");
1469  NumberFree(Num3,"modinverses");
1470  NumberFree(Num1,"modinverses");
1471  NumberFree(Num2,"modinverses");
1472  t[2] = 0; /* mark function as clean. */
1473  goto Restart;
1474  }
1475  break;
1476  case GCDFUNCTION:
1477 #ifdef EVALUATEGCD
1478 #ifdef NEWGCDFUNCTION
1479  {
1480 /*
1481  Has two integer arguments
1482  Four cases: S,S, S,L, L,S, L,L
1483 */
1484  WORD *num1, *num2, size1, size2, stor1, stor2, *ttt, ti;
1485  if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1486  && t[FUNHEAD+2] == -SNUMBER && t[FUNHEAD+1] != 0
1487  && t[FUNHEAD+3] != 0 ) { /* Short,Short */
1488  stor1 = t[FUNHEAD+1];
1489  stor2 = t[FUNHEAD+3];
1490  if ( stor1 < 0 ) stor1 = -stor1;
1491  if ( stor2 < 0 ) stor2 = -stor2;
1492  num1 = &stor1; num2 = &stor2;
1493  size1 = size2 = 1;
1494  goto gcdcalc;
1495  }
1496  else if ( t[1] > FUNHEAD+4 ) {
1497  if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] != 0
1498  && t[FUNHEAD+2] == t[1]-FUNHEAD-2 &&
1499  ABS(t[t[1]-1]) == t[FUNHEAD+2]-1-ARGHEAD ) { /* Short,Long */
1500  num2 = t + t[1];
1501  size2 = ABS(num2[-1]);
1502  ttt = num2-1;
1503  num2 -= size2;
1504  size2 = (size2-1)/2;
1505  ti = size2;
1506  while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1507  if ( ti == 1 && ttt[-1] == 1 ) {
1508  stor1 = t[FUNHEAD+1];
1509  if ( stor1 < 0 ) stor1 = -stor1;
1510  num1 = &stor1;
1511  size1 = 1;
1512  goto gcdcalc;
1513  }
1514  else pcom[ncom++] = t;
1515  }
1516  else if ( t[FUNHEAD] > 0 &&
1517  t[FUNHEAD]-1-ARGHEAD == ABS(t[t[FUNHEAD]+FUNHEAD-1]) ) {
1518  num1 = t + FUNHEAD + t[FUNHEAD];
1519  size1 = ABS(num1[-1]);
1520  ttt = num1-1;
1521  num1 -= size1;
1522  size1 = (size1-1)/2;
1523  ti = size1;
1524  while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1525  if ( ti == 1 && ttt[-1] == 1 ) {
1526  if ( t[1]-FUNHEAD == t[FUNHEAD]+2 && t[t[1]-2] == -SNUMBER
1527  && t[t[1]-1] != 0 ) { /* Long,Short */
1528  stor2 = t[t[1]-1];
1529  if ( stor2 < 0 ) stor2 = -stor2;
1530  num2 = &stor2;
1531  size2 = 1;
1532  goto gcdcalc;
1533  }
1534  else if ( t[1]-FUNHEAD == t[FUNHEAD]+t[FUNHEAD+t[FUNHEAD]]
1535  && ABS(t[t[1]-1]) == t[FUNHEAD+t[FUNHEAD]] - ARGHEAD-1 ) {
1536  num2 = t + t[1];
1537  size2 = ABS(num2[-1]);
1538  ttt = num2-1;
1539  num2 -= size2;
1540  size2 = (size2-1)/2;
1541  ti = size2;
1542  while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1543  if ( ti == 1 && ttt[-1] == 1 ) {
1544 gcdcalc: if ( GcdLong(BHEAD (UWORD *)num1,size1,(UWORD *)num2,size2
1545  ,(UWORD *)lnum,&nnum) ) goto FromNorm;
1546  goto MulIn;
1547  }
1548  else pcom[ncom++] = t;
1549  }
1550  else pcom[ncom++] = t;
1551  }
1552  else pcom[ncom++] = t;
1553  }
1554  else pcom[ncom++] = t;
1555  }
1556  else pcom[ncom++] = t;
1557  }
1558 #else
1559  {
1560  WORD *gcd = AT.WorkPointer;
1561  if ( ( gcd = EvaluateGcd(BHEAD t) ) == 0 ) goto FromNorm;
1562  if ( *gcd == 4 && gcd[1] == 1 && gcd[2] == 1 && gcd[4] == 0 ) {
1563  AT.WorkPointer = gcd;
1564  }
1565  else if ( gcd[*gcd] == 0 ) {
1566  WORD *t1, iii, change, *num, *den, numsize, densize;
1567  if ( gcd[*gcd-1] < *gcd-1 ) {
1568  t1 = gcd+1;
1569  for ( iii = 2; iii < t1[1]; iii += 2 ) {
1570  change = ExtraSymbol(t1[iii],t1[iii+1],nsym,ppsym,&ncoef);
1571  nsym += change;
1572  ppsym += change << 1;
1573  }
1574  }
1575  t1 = gcd + *gcd;
1576  iii = t1[-1]; num = t1-iii; numsize = (iii-1)/2;
1577  den = num + numsize; densize = numsize;
1578  while ( numsize > 1 && num[numsize-1] == 0 ) numsize--;
1579  while ( densize > 1 && den[densize-1] == 0 ) densize--;
1580  if ( numsize > 1 || num[0] != 1 ) {
1581  ncoef = REDLENG(ncoef);
1582  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)num,numsize) ) goto FromNorm;
1583  ncoef = INCLENG(ncoef);
1584  }
1585  if ( densize > 1 || den[0] != 1 ) {
1586  ncoef = REDLENG(ncoef);
1587  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)den,densize) ) goto FromNorm;
1588  ncoef = INCLENG(ncoef);
1589  }
1590  AT.WorkPointer = gcd;
1591  }
1592  else { /* a whole expression */
1593 /*
1594  Action: Put the expression in a compiler buffer.
1595  Insert a SUBEXPRESSION subterm
1596  Set the return value of the routine such that in
1597  Generator the term gets sent again to TestSub.
1598 
1599  1: put in C (ebufnum)
1600  2: after that the WorkSpace is free again.
1601  3: insert the SUBEXPRESSION
1602  4: copy the top part of the term down
1603 */
1604  LONG size = AT.WorkPointer - gcd;
1605 
1606  ss = AddRHS(AT.ebufnum,1);
1607  while ( (ss + size + 10) > C->Top ) ss = DoubleCbuffer(AT.ebufnum,ss);
1608  tt = gcd;
1609  NCOPY(ss,tt,size);
1610  C->rhs[C->numrhs+1] = ss;
1611  C->Pointer = ss;
1612 
1613  t[0] = SUBEXPRESSION;
1614  t[1] = SUBEXPSIZE;
1615  t[2] = C->numrhs;
1616  t[3] = 1;
1617  t[4] = AT.ebufnum;
1618  t += 5;
1619  tt = term + *term;
1620  while ( r < tt ) *t++ = *r++;
1621  *term = t - term;
1622 
1623  regval = 1;
1624  goto RegEnd;
1625  }
1626  }
1627 #endif
1628 #else
1629  MesPrint(" Unexpected call to EvaluateGCD");
1630  Terminate(-1);
1631 #endif
1632  break;
1633  case MINFUNCTION:
1634  case MAXFUNCTION:
1635  if ( t[1] == FUNHEAD ) break;
1636  {
1637  WORD *ttt = t + FUNHEAD;
1638  WORD *tttstop = t + t[1];
1639  WORD tterm[4], iii;
1640  while ( ttt < tttstop ) {
1641  if ( *ttt > 0 ) {
1642  if ( ttt[ARGHEAD]-1 > ABS(ttt[*ttt-1]) ) goto nospec;
1643  ttt += *ttt;
1644  }
1645  else {
1646  if ( *ttt != -SNUMBER ) goto nospec;
1647  ttt += 2;
1648  }
1649  }
1650 /*
1651  Function has only numerical arguments
1652  Pick up the first argument.
1653 */
1654  ttt = t + FUNHEAD;
1655  if ( *ttt > 0 ) {
1656 loadnew1:
1657  for ( iii = 0; iii < ttt[ARGHEAD]; iii++ )
1658  AT.n_llnum[iii] = ttt[ARGHEAD+iii];
1659  ttt += *ttt;
1660  }
1661  else {
1662 loadnew2:
1663  if ( ttt[1] == 0 ) {
1664  AT.n_llnum[0] = AT.n_llnum[1] = AT.n_llnum[2] = AT.n_llnum[3] = 0;
1665  }
1666  else {
1667  AT.n_llnum[0] = 4;
1668  if ( ttt[1] > 0 ) { AT.n_llnum[1] = ttt[1]; AT.n_llnum[3] = 3; }
1669  else { AT.n_llnum[1] = -ttt[1]; AT.n_llnum[3] = -3; }
1670  AT.n_llnum[2] = 1;
1671  }
1672  ttt += 2;
1673  }
1674 /*
1675  Now loop over the other arguments
1676 */
1677  while ( ttt < tttstop ) {
1678  if ( *ttt > 0 ) {
1679  if ( AT.n_llnum[0] == 0 ) {
1680  if ( ( *t == MINFUNCTION && ttt[*ttt-1] < 0 )
1681  || ( *t == MAXFUNCTION && ttt[*ttt-1] > 0 ) )
1682  goto loadnew1;
1683  }
1684  else {
1685  ttt += ARGHEAD;
1686  iii = CompCoef(AT.n_llnum,ttt);
1687  if ( ( iii > 0 && *t == MINFUNCTION )
1688  || ( iii < 0 && *t == MAXFUNCTION ) ) {
1689  for ( iii = 0; iii < ttt[0]; iii++ )
1690  AT.n_llnum[iii] = ttt[iii];
1691  }
1692  }
1693  ttt += *ttt;
1694  }
1695  else {
1696  if ( AT.n_llnum[0] == 0 ) {
1697  if ( ( *t == MINFUNCTION && ttt[1] < 0 )
1698  || ( *t == MAXFUNCTION && ttt[1] > 0 ) )
1699  goto loadnew2;
1700  }
1701  else if ( ttt[1] == 0 ) {
1702  if ( ( *t == MINFUNCTION && AT.n_llnum[*AT.n_llnum-1] > 0 )
1703  || ( *t == MAXFUNCTION && AT.n_llnum[*AT.n_llnum-1] < 0 ) ) {
1704  AT.n_llnum[0] = 0;
1705  }
1706  }
1707  else {
1708  tterm[0] = 4; tterm[2] = 1;
1709  if ( ttt[1] < 0 ) { tterm[1] = -ttt[1]; tterm[3] = -3; }
1710  else { tterm[1] = ttt[1]; tterm[3] = 3; }
1711  iii = CompCoef(AT.n_llnum,tterm);
1712  if ( ( iii > 0 && *t == MINFUNCTION )
1713  || ( iii < 0 && *t == MAXFUNCTION ) ) {
1714  for ( iii = 0; iii < 4; iii++ )
1715  AT.n_llnum[iii] = tterm[iii];
1716  }
1717  }
1718  ttt += 2;
1719  }
1720  }
1721  if ( AT.n_llnum[0] == 0 ) goto NormZero;
1722  ncoef = REDLENG(ncoef);
1723  nnum = REDLENG(AT.n_llnum[*AT.n_llnum-1]);
1724  if ( MulRat(BHEAD (UWORD *)AT.n_coef,ncoef,(UWORD *)lnum,nnum,
1725  (UWORD *)AT.n_coef,&ncoef) ) goto FromNorm;
1726  ncoef = INCLENG(ncoef);
1727  }
1728  break;
1729  case INVERSEFACTORIAL:
1730  if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] >= 0 ) {
1731  if ( Factorial(BHEAD t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
1732  goto FromNorm;
1733  ncoef = REDLENG(ncoef);
1734  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) ) goto FromNorm;
1735  ncoef = INCLENG(ncoef);
1736  }
1737  else {
1738 nospec: pcom[ncom++] = t;
1739  }
1740  break;
1741  case MAXPOWEROF:
1742  if ( ( t[FUNHEAD] == -SYMBOL )
1743  && ( t[FUNHEAD+1] > 0 ) && ( t[1] == FUNHEAD+2 ) ) {
1744  *((UWORD *)lnum) = symbols[t[FUNHEAD+1]].maxpower;
1745  nnum = 1;
1746  goto MulIn;
1747  }
1748  else { pcom[ncom++] = t; }
1749  break;
1750  case MINPOWEROF:
1751  if ( ( t[FUNHEAD] == -SYMBOL )
1752  && ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+2 ) ) {
1753  *((UWORD *)lnum) = symbols[t[FUNHEAD+1]].minpower;
1754  nnum = 1;
1755  goto MulIn;
1756  }
1757  else { pcom[ncom++] = t; }
1758  break;
1759  case PRIMENUMBER :
1760  if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER
1761  && t[FUNHEAD+1] > 0 ) {
1762  UWORD xp = (UWORD)(NextPrime(BHEAD t[FUNHEAD+1]));
1763  ncoef = REDLENG(ncoef);
1764  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,&xp,1) ) goto FromNorm;
1765  ncoef = INCLENG(ncoef);
1766  }
1767  else goto defaultcase;
1768  break;
1769  case LNUMBER :
1770  ncoef = REDLENG(ncoef);
1771  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)(t+3),t[2]) ) goto FromNorm;
1772  ncoef = INCLENG(ncoef);
1773  break;
1774  case SNUMBER :
1775  if ( t[2] < 0 ) {
1776  t[2] = -t[2];
1777  if ( t[3] & 1 ) ncoef = -ncoef;
1778  }
1779  else if ( t[2] == 0 ) {
1780  if ( t[3] < 0 ) goto NormInf;
1781  goto NormZero;
1782  }
1783  lnum[0] = t[2];
1784  nnum = 1;
1785  if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) ) goto FromNorm;
1786  ncoef = REDLENG(ncoef);
1787  if ( t[3] < 0 ) {
1788  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) )
1789  goto FromNorm;
1790  }
1791  else if ( t[3] > 0 ) {
1792  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) )
1793  goto FromNorm;
1794  }
1795  ncoef = INCLENG(ncoef);
1796  break;
1797  case GAMMA :
1798  case GAMMAI :
1799  case GAMMAFIVE :
1800  case GAMMASIX :
1801  case GAMMASEVEN :
1802  pnco[nnco++] = t;
1803  t += FUNHEAD+1;
1804  goto ScanCont;
1805  case LEVICIVITA :
1806  peps[neps++] = t;
1807  if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG ) {
1808  t[2] &= ~DIRTYFLAG;
1809  t[2] |= DIRTYSYMFLAG;
1810  }
1811  t += FUNHEAD;
1812 ScanCont: while ( t < r ) {
1813  if ( *t >= AM.OffsetIndex &&
1814  ( *t >= AM.DumInd || ( *t < AM.WilInd &&
1815  indices[*t-AM.OffsetIndex].dimension ) ) )
1816  pcon[ncon++] = t;
1817  t++;
1818  }
1819  break;
1820  case EXPONENT :
1821  {
1822  WORD *rr;
1823  k = 1;
1824  rr = t + FUNHEAD;
1825  if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) )
1826  k = 0;
1827  if ( *rr == -SNUMBER && rr[1] == 1 ) break;
1828  if ( *rr <= -FUNCTION ) k = *rr;
1829  NEXTARG(rr)
1830  if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) ) {
1831  if ( k == 0 ) goto NormZZ;
1832  break;
1833  }
1834  if ( *rr == -SNUMBER && rr[1] > 0 && rr[1] < MAXPOWER && k < 0 ) {
1835  k = -k;
1836  if ( functions[k-FUNCTION].commute ) {
1837  for ( i = 0; i < rr[1]; i++ ) pnco[nnco++] = rr-1;
1838  }
1839  else {
1840  for ( i = 0; i < rr[1]; i++ ) pcom[ncom++] = rr-1;
1841  }
1842  break;
1843  }
1844  if ( k == 0 ) goto NormZero;
1845  if ( t[FUNHEAD] == -SYMBOL && *rr == -SNUMBER && t[1] == FUNHEAD+4 ) {
1846  if ( rr[1] < MAXPOWER ) {
1847  t = rr; *rr = t[FUNHEAD+1];
1848  goto NextSymbol;
1849  }
1850  }
1851 /*
1852  if ( ( t[FUNHEAD] > 0 && t[FUNHEAD+1] != 0 )
1853  || ( *rr > 0 && rr[1] != 0 ) ) {}
1854  else
1855 */
1856  t[2] &= ~DIRTYSYMFLAG;
1857 
1858  pnco[nnco++] = t;
1859  }
1860  break;
1861  case DENOMINATOR :
1862  t[2] &= ~DIRTYSYMFLAG;
1863  pden[nden++] = t;
1864  pnco[nnco++] = t;
1865  break;
1866  case INDEX :
1867  t += 2;
1868  do {
1869  if ( *t == 0 ) goto NormZero;
1870  if ( *t > 0 && *t < AM.OffsetIndex ) {
1871  lnum[0] = *t++;
1872  nnum = 1;
1873  ncoef = REDLENG(ncoef);
1874  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)lnum,nnum) )
1875  goto FromNorm;
1876  ncoef = INCLENG(ncoef);
1877  }
1878  else if ( *t == NOINDEX ) t++;
1879  else pind[nind++] = *t++;
1880  } while ( t < r );
1881  break;
1882  case SUBEXPRESSION :
1883  if ( t[3] == 0 ) break;
1884  case EXPRESSION :
1885  goto RegEnd;
1886  case ROOTFUNCTION :
1887 /*
1888  Tries to take the n-th root inside the rationals
1889  If this is not possible, it clears all flags and
1890  hence tries no more.
1891  Notation:
1892  root_(power(=integer),(rational)number)
1893 */
1894  { WORD nc;
1895  if ( t[2] == 0 ) goto defaultcase;
1896  if ( t[FUNHEAD] != -SNUMBER || t[FUNHEAD+1] < 0 ) goto defaultcase;
1897  if ( t[FUNHEAD+2] == -SNUMBER ) {
1898  if ( t[FUNHEAD+1] == 0 && t[FUNHEAD+3] == 0 ) goto NormZZ;
1899  if ( t[FUNHEAD+1] == 0 ) break;
1900  if ( t[FUNHEAD+3] < 0 ) {
1901  AT.WorkPointer[0] = -t[FUNHEAD+3];
1902  nc = -1;
1903  }
1904  else {
1905  AT.WorkPointer[0] = t[FUNHEAD+3];
1906  nc = 1;
1907  }
1908  AT.WorkPointer[1] = 1;
1909  }
1910  else if ( t[FUNHEAD+2] == t[1]-FUNHEAD-2
1911  && t[FUNHEAD+2] == t[FUNHEAD+2+ARGHEAD]+ARGHEAD
1912  && ABS(t[t[1]-1]) == t[FUNHEAD+2+ARGHEAD] - 1 ) {
1913  WORD *r1, *r2;
1914  if ( t[FUNHEAD+1] == 0 ) break;
1915  i = t[t[1]-1]; r1 = t + FUNHEAD+ARGHEAD+3;
1916  nc = REDLENG(i);
1917  i = ABS(i) - 1;
1918  r2 = AT.WorkPointer;
1919  while ( --i >= 0 ) *r2++ = *r1++;
1920  }
1921  else goto defaultcase;
1922  if ( TakeRatRoot((UWORD *)AT.WorkPointer,&nc,t[FUNHEAD+1]) ) {
1923  t[2] = 0;
1924  goto defaultcase;
1925  }
1926  ncoef = REDLENG(ncoef);
1927  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)AT.WorkPointer,nc) )
1928  goto FromNorm;
1929  if ( nc < 0 ) nc = -nc;
1930  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)(AT.WorkPointer+nc),nc) )
1931  goto FromNorm;
1932  ncoef = INCLENG(ncoef);
1933  }
1934  break;
1935  case RANDOMFUNCTION :
1936  {
1937  WORD nnc, nc, nca, nr;
1938  UWORD xx;
1939 /*
1940  Needs one positive integer argument.
1941  returns (wranf()%argument)+1.
1942  We may call wranf several times to paste UWORDS together
1943  when we need long numbers.
1944  We make little errors when taking the % operator
1945  (not 100% uniform). We correct for that by redoing the
1946  calculation in the (unlikely) case that we are in leftover area
1947 */
1948  if ( t[1] == FUNHEAD ) goto defaultcase;
1949  if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER &&
1950  t[FUNHEAD+1] > 0 ) {
1951  if ( t[FUNHEAD+1] == 1 ) break;
1952 redoshort:
1953  ((UWORD *)AT.WorkPointer)[0] = wranf(BHEAD0);
1954  ((UWORD *)AT.WorkPointer)[1] = wranf(BHEAD0);
1955  nr = 2;
1956  if ( ((UWORD *)AT.WorkPointer)[1] == 0 ) {
1957  nr = 1;
1958  if ( ((UWORD *)AT.WorkPointer)[0] == 0 ) {
1959  nr = 0;
1960  }
1961  }
1962  xx = (UWORD)(t[FUNHEAD+1]);
1963  if ( nr ) {
1964  DivLong((UWORD *)AT.WorkPointer,nr
1965  ,&xx,1
1966  ,((UWORD *)AT.WorkPointer)+4,&nnc
1967  ,((UWORD *)AT.WorkPointer)+2,&nc);
1968  ((UWORD *)AT.WorkPointer)[4] = 0;
1969  ((UWORD *)AT.WorkPointer)[5] = 0;
1970  ((UWORD *)AT.WorkPointer)[6] = 1;
1971  DivLong((UWORD *)AT.WorkPointer+4,3
1972  ,&xx,1
1973  ,((UWORD *)AT.WorkPointer)+9,&nnc
1974  ,((UWORD *)AT.WorkPointer)+7,&nca);
1975  AddLong((UWORD *)AT.WorkPointer+4,3
1976  ,((UWORD *)AT.WorkPointer)+7,-nca
1977  ,((UWORD *)AT.WorkPointer)+9,&nnc);
1978  if ( BigLong((UWORD *)AT.WorkPointer,nr
1979  ,((UWORD *)AT.WorkPointer)+9,nnc) >= 0 ) goto redoshort;
1980  }
1981  else nc = 0;
1982  if ( nc == 0 ) {
1983  AT.WorkPointer[2] = (WORD)xx;
1984  nc = 1;
1985  }
1986  ncoef = REDLENG(ncoef);
1987  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+2,nc) )
1988  goto FromNorm;
1989  ncoef = INCLENG(ncoef);
1990  }
1991  else if ( t[FUNHEAD] > 0 && t[1] == t[FUNHEAD]+FUNHEAD
1992  && ABS(t[t[1]-1]) == t[FUNHEAD]-1-ARGHEAD && t[t[1]-1] > 0 ) {
1993  WORD nna, nnb, nni, nnb2, nnb2a;
1994  UWORD *nnt;
1995  nna = t[t[1]-1];
1996  nnb2 = nna-1;
1997  nnb = nnb2/2;
1998  nnt = (UWORD *)(t+t[1]-1-nnb); /* start of denominator */
1999  if ( *nnt != 1 ) goto defaultcase;
2000  for ( nni = 1; nni < nnb; nni++ ) {
2001  if ( nnt[nni] != 0 ) goto defaultcase;
2002  }
2003  nnt = (UWORD *)(t + FUNHEAD + ARGHEAD + 1);
2004 
2005  for ( nni = 0; nni < nnb2; nni++ ) {
2006  ((UWORD *)AT.WorkPointer)[nni] = wranf(BHEAD0);
2007  }
2008  nnb2a = nnb2;
2009  while ( nnb2a > 0 && ((UWORD *)AT.WorkPointer)[nnb2a-1] == 0 ) nnb2a--;
2010  if ( nnb2a > 0 ) {
2011  DivLong((UWORD *)AT.WorkPointer,nnb2a
2012  ,nnt,nnb
2013  ,((UWORD *)AT.WorkPointer)+2*nnb2,&nnc
2014  ,((UWORD *)AT.WorkPointer)+nnb2,&nc);
2015  for ( nni = 0; nni < nnb2; nni++ ) {
2016  ((UWORD *)AT.WorkPointer)[nni+2*nnb2] = 0;
2017  }
2018  ((UWORD *)AT.WorkPointer)[3*nnb2] = 1;
2019  DivLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2020  ,nnt,nnb
2021  ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc
2022  ,((UWORD *)AT.WorkPointer)+3*nnb2+1,&nca);
2023  AddLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2024  ,((UWORD *)AT.WorkPointer)+3*nnb2+1,-nca
2025  ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc);
2026  if ( BigLong((UWORD *)AT.WorkPointer,nnb2a
2027  ,((UWORD *)AT.WorkPointer)+4*nnb2+1,nnc) >= 0 ) goto redoshort;
2028  }
2029  else nc = 0;
2030  if ( nc == 0 ) {
2031  for ( nni = 0; nni < nnb; nni++ ) {
2032  ((UWORD *)AT.WorkPointer)[nnb2+nni] = nnt[nni];
2033  }
2034  nc = nnb;
2035  }
2036  ncoef = REDLENG(ncoef);
2037  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+nnb2,nc) )
2038  goto FromNorm;
2039  ncoef = INCLENG(ncoef);
2040  }
2041  else goto defaultcase;
2042  }
2043  break;
2044  case RANPERM :
2045  if ( *t == RANPERM && t[1] > FUNHEAD && t[FUNHEAD] <= -FUNCTION ) {
2046  WORD **pwork;
2047  WORD *mm, *ww, *ow = AT.WorkPointer;
2048  WORD *Array, *targ, *argstop, narg = 0, itot;
2049  int ie;
2050  argstop = t+t[1];
2051  targ = t+FUNHEAD+1;
2052  while ( targ < argstop ) {
2053  narg++; NEXTARG(targ);
2054  }
2055  WantAddPointers(narg);
2056  pwork = AT.pWorkSpace + AT.pWorkPointer;
2057  targ = t+FUNHEAD+1; narg = 0;
2058  while ( targ < argstop ) {
2059  pwork[narg++] = targ;
2060  NEXTARG(targ);
2061  }
2062 /*
2063  Make a random permutation of the numbers 0,...,narg-1
2064  The following code works also for narg == 0 and narg == 1
2065 */
2066  ow = AT.WorkPointer;
2067  Array = AT.WorkPointer;
2068  AT.WorkPointer += narg;
2069  for ( i = 0; i < narg; i++ ) Array[i] = i;
2070  for ( i = 2; i <= narg; i++ ) {
2071  itot = (WORD)(iranf(BHEAD i));
2072  for ( j = 0; j < itot; j++ ) CYCLE1(WORD,Array,i)
2073  }
2074  mm = AT.WorkPointer;
2075  *mm++ = -t[FUNHEAD];
2076  *mm++ = t[1] - 1;
2077  for ( ie = 2; ie < FUNHEAD; ie++ ) *mm++ = t[ie];
2078  for ( i = 0; i < narg; i++ ) {
2079  ww = pwork[Array[i]];
2080  CopyArg(mm,ww);
2081  }
2082  mm = AT.WorkPointer; t++; ww = t;
2083  i = mm[1]; NCOPY(ww,mm,i)
2084  AT.WorkPointer = ow;
2085  goto TryAgain;
2086  }
2087  pnco[nnco++] = t;
2088  break;
2089  case INTFUNCTION :
2090 /*
2091  Can be resolved if the first argument is a number
2092  and the second argument either doesn't exist or has
2093  the value +1, 0, -1
2094  +1 : rounding up
2095  0 : rounding towards zero
2096  -1 : rounding down (same as no argument)
2097 */
2098  if ( t[1] <= FUNHEAD ) break;
2099  {
2100  WORD *rr, den, num;
2101  to = t + FUNHEAD;
2102  if ( *to > 0 ) {
2103  if ( *to == ARGHEAD ) goto NormZero;
2104  rr = to + *to;
2105  i = rr[-1];
2106  j = ABS(i);
2107  if ( to[ARGHEAD] != j+1 ) goto NoInteg;
2108  if ( rr >= r ) k = -1;
2109  else if ( *rr == ARGHEAD ) { k = 0; rr += ARGHEAD; }
2110  else if ( *rr == -SNUMBER ) { k = rr[1]; rr += 2; }
2111  else goto NoInteg;
2112  if ( rr != r ) goto NoInteg;
2113  if ( k > 1 || k < -1 ) goto NoInteg;
2114  to += ARGHEAD+1;
2115  j = (j-1) >> 1;
2116  i = ( i < 0 ) ? -j: j;
2117  UnPack((UWORD *)to,i,&den,&num);
2118 /*
2119  Potentially the use of NoScrat2 is unsafe.
2120  It makes the routine not reentrant, but because it is
2121  used only locally and because we only call the
2122  low level routines DivLong and AddLong which never
2123  make calls involving Normalize, things are OK after all
2124 */
2125  if ( AN.NoScrat2 == 0 ) {
2126  AN.NoScrat2 = (UWORD *)Malloc1((AM.MaxTal+2)*sizeof(UWORD),"Normalize");
2127  }
2128  if ( DivLong((UWORD *)to,num,(UWORD *)(to+j),den
2129  ,(UWORD *)AT.WorkPointer,&num,AN.NoScrat2,&den) ) goto FromNorm;
2130  if ( k < 0 && den < 0 ) {
2131  *AN.NoScrat2 = 1;
2132  den = -1;
2133  if ( AddLong((UWORD *)AT.WorkPointer,num
2134  ,AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2135  goto FromNorm;
2136  }
2137  else if ( k > 0 && den > 0 ) {
2138  *AN.NoScrat2 = 1;
2139  den = 1;
2140  if ( AddLong((UWORD *)AT.WorkPointer,num,
2141  AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2142  goto FromNorm;
2143  }
2144 
2145  }
2146  else if ( *to == -SNUMBER ) { /* No rounding needed */
2147  if ( to[1] < 0 ) { *AT.WorkPointer = -to[1]; num = -1; }
2148  else if ( to[1] == 0 ) goto NormZero;
2149  else { *AT.WorkPointer = to[1]; num = 1; }
2150  }
2151  else goto NoInteg;
2152  if ( num == 0 ) goto NormZero;
2153  ncoef = REDLENG(ncoef);
2154  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)AT.WorkPointer,num) )
2155  goto FromNorm;
2156  ncoef = INCLENG(ncoef);
2157  break;
2158  }
2159 NoInteg:;
2160 /*
2161  Fall through if it cannot be resolved
2162 */
2163  default :
2164 defaultcase:;
2165  if ( *t < FUNCTION ) {
2166  MLOCK(ErrorMessageLock);
2167  MesPrint("Illegal code in Norm");
2168 #ifdef DEBUGON
2169  {
2170  UBYTE OutBuf[140];
2171  AO.OutFill = AO.OutputLine = OutBuf;
2172  t = term;
2173  AO.OutSkip = 3;
2174  FiniLine();
2175  i = *t;
2176  while ( --i >= 0 ) {
2177  TalToLine((UWORD)(*t++));
2178  TokenToLine((UBYTE *)" ");
2179  }
2180  AO.OutSkip = 0;
2181  FiniLine();
2182  }
2183 #endif
2184  MUNLOCK(ErrorMessageLock);
2185  goto NormMin;
2186  }
2187  if ( *t == REPLACEMENT && ReplaceType == -1 ) {
2188  ReplaceType = 0;
2189  from = t;
2190  if ( AN.RSsize < 2*t[1]+SUBEXPSIZE ) {
2191  if ( AN.ReplaceScrat ) M_free(AN.ReplaceScrat,"AN.ReplaceScrat");
2192  AN.RSsize = 2*t[1]+SUBEXPSIZE+40;
2193  AN.ReplaceScrat = (WORD *)Malloc1((AN.RSsize+1)*sizeof(WORD),"AN.ReplaceScrat");
2194  }
2195  t += FUNHEAD;
2196  ReplaceSub = AN.ReplaceScrat;
2197  ReplaceSub += SUBEXPSIZE;
2198  while ( t < r ) {
2199  if ( *t > 0 ) goto NoRep;
2200  if ( *t <= -FUNCTION ) {
2201  *ReplaceSub++ = FUNTOFUN;
2202  *ReplaceSub++ = 4;
2203  *ReplaceSub++ = -*t++;
2204  if ( *t > -FUNCTION ) goto NoRep;
2205  *ReplaceSub++ = -*t++;
2206  }
2207  else if ( t+4 > r ) goto NoRep;
2208  else {
2209  if ( *t == -SYMBOL ) {
2210  if ( t[2] == -SYMBOL && t+4 <= r )
2211  *ReplaceSub++ = SYMTOSYM;
2212  else if ( t[2] == -SNUMBER && t+4 <= r ) {
2213  *ReplaceSub++ = SYMTONUM;
2214  if ( ReplaceType == 0 ) {
2215  oldtoprhs = C->numrhs;
2216  oldcpointer = C->Pointer - C->Buffer;
2217  }
2218  ReplaceType = 1;
2219  }
2220  else if ( t[2] == ARGHEAD && t+2+ARGHEAD <= r ) {
2221  *ReplaceSub++ = SYMTONUM;
2222  *ReplaceSub++ = 4;
2223  *ReplaceSub++ = t[1];
2224  *ReplaceSub++ = 0;
2225  t += 2+ARGHEAD;
2226  continue;
2227  }
2228 /*
2229  Next is the subexpression. We have to test that
2230  it isn't vector-like or index-like
2231 */
2232  else if ( t[2] > 0 ) {
2233  WORD *sstop, *ttstop, n;
2234  ss = t+2;
2235  sstop = ss + *ss;
2236  ss += ARGHEAD;
2237  while ( ss < sstop ) {
2238  tt = ss + *ss;
2239  ttstop = tt - ABS(tt[-1]);
2240  ss++;
2241  while ( ss < ttstop ) {
2242  if ( *ss == INDEX ) goto NoRep;
2243  ss += ss[1];
2244  }
2245  ss = tt;
2246  }
2247  subtype = SYMTOSUB;
2248  if ( ReplaceType == 0 ) {
2249  oldtoprhs = C->numrhs;
2250  oldcpointer = C->Pointer - C->Buffer;
2251  }
2252  ReplaceType = 1;
2253  ss = AddRHS(AT.ebufnum,1);
2254  tt = t+2;
2255  n = *tt - ARGHEAD;
2256  tt += ARGHEAD;
2257  while ( (ss + n + 10) > C->Top ) ss = DoubleCbuffer(AT.ebufnum,ss);
2258  while ( --n >= 0 ) *ss++ = *tt++;
2259  *ss++ = 0;
2260  C->rhs[C->numrhs+1] = ss;
2261  C->Pointer = ss;
2262  *ReplaceSub++ = subtype;
2263  *ReplaceSub++ = 4;
2264  *ReplaceSub++ = t[1];
2265  *ReplaceSub++ = C->numrhs;
2266  t += 2 + t[2];
2267  continue;
2268  }
2269  else goto NoRep;
2270  }
2271  else if ( *t == -VECTOR && t+4 <= r ) {
2272  if ( t[2] == -VECTOR ) *ReplaceSub++ = VECTOVEC;
2273  else if ( t[2] == -MINVECTOR )
2274  *ReplaceSub++ = VECTOMIN;
2275 /*
2276  Next is a vector-like subexpression
2277  Search for vector nature first
2278 */
2279  else if ( t[2] > 0 ) {
2280  WORD *sstop, *ttstop, *w, *mm, n, count;
2281  WORD *v1, *v2 = 0;
2282  ss = t+2;
2283  sstop = ss + *ss;
2284  ss += ARGHEAD;
2285  while ( ss < sstop ) {
2286  tt = ss + *ss;
2287  ttstop = tt - ABS(tt[-1]);
2288  ss++;
2289  count = 0;
2290  while ( ss < ttstop ) {
2291  if ( *ss == INDEX ) {
2292  n = ss[1] - 2; ss += 2;
2293  while ( --n >= 0 ) {
2294  if ( *ss < MINSPEC ) count++;
2295  ss++;
2296  }
2297  }
2298  else ss += ss[1];
2299  }
2300  if ( count != 1 ) goto NoRep;
2301  ss = tt;
2302  }
2303  subtype = VECTOSUB;
2304  if ( ReplaceType == 0 ) {
2305  oldtoprhs = C->numrhs;
2306  oldcpointer = C->Pointer - C->Buffer;
2307  }
2308  ReplaceType = 1;
2309  mm = AddRHS(AT.ebufnum,1);
2310  *ReplaceSub++ = subtype;
2311  *ReplaceSub++ = 4;
2312  *ReplaceSub++ = t[1];
2313  *ReplaceSub++ = C->numrhs;
2314  w = t+2;
2315  n = *w - ARGHEAD;
2316  w += ARGHEAD;
2317  while ( (mm + n + 10) > C->Top )
2318  mm = DoubleCbuffer(AT.ebufnum,mm);
2319  while ( --n >= 0 ) *mm++ = *w++;
2320  *mm++ = 0;
2321  C->rhs[C->numrhs+1] = mm;
2322  C->Pointer = mm;
2323  mm = AddRHS(AT.ebufnum,1);
2324  w = t+2;
2325  n = *w - ARGHEAD;
2326  w += ARGHEAD;
2327  while ( (mm + n + 13) > C->Top )
2328  mm = DoubleCbuffer(AT.ebufnum,mm);
2329  sstop = w + n;
2330  while ( w < sstop ) {
2331  tt = w + *w; ttstop = tt - ABS(tt[-1]);
2332  ss = mm; mm++; w++;
2333  while ( w < ttstop ) { /* Subterms */
2334  if ( *w != INDEX ) {
2335  n = w[1];
2336  NCOPY(mm,w,n);
2337  }
2338  else {
2339  v1 = mm;
2340  *mm++ = *w++;
2341  *mm++ = n = *w++;
2342  n -= 2;
2343  while ( --n >= 0 ) {
2344  if ( *w >= MINSPEC ) *mm++ = *w++;
2345  else v2 = w++;
2346  }
2347  n = WORDDIF(mm,v1);
2348  if ( n != v1[1] ) {
2349  if ( n <= 2 ) mm -= 2;
2350  else v1[1] = n;
2351  *mm++ = VECTOR;
2352  *mm++ = 4;
2353  *mm++ = *v2;
2354  *mm++ = FUNNYVEC;
2355  }
2356  }
2357  }
2358  while ( w < tt ) *mm++ = *w++;
2359  *ss = WORDDIF(mm,ss);
2360  }
2361  *mm++ = 0;
2362  C->rhs[C->numrhs+1] = mm;
2363  C->Pointer = mm;
2364  if ( mm > C->Top ) {
2365  MLOCK(ErrorMessageLock);
2366  MesPrint("Internal error in Normalize with extra compiler buffer");
2367  MUNLOCK(ErrorMessageLock);
2368  Terminate(-1);
2369  }
2370  t += 2 + t[2];
2371  continue;
2372  }
2373  else goto NoRep;
2374  }
2375  else if ( *t == -INDEX ) {
2376  if ( ( t[2] == -INDEX || t[2] == -VECTOR )
2377  && t+4 <= r )
2378  *ReplaceSub++ = INDTOIND;
2379  else if ( t[1] >= AM.OffsetIndex ) {
2380  if ( t[2] == -SNUMBER && t+4 <= r
2381  && t[3] >= 0 && t[3] < AM.OffsetIndex )
2382  *ReplaceSub++ = INDTOIND;
2383  else if ( t[2] == ARGHEAD && t+2+ARGHEAD <= r ) {
2384  *ReplaceSub++ = INDTOIND;
2385  *ReplaceSub++ = 4;
2386  *ReplaceSub++ = t[1];
2387  *ReplaceSub++ = 0;
2388  t += 2+ARGHEAD;
2389  continue;
2390  }
2391  else goto NoRep;
2392  }
2393  else goto NoRep;
2394  }
2395  else goto NoRep;
2396  *ReplaceSub++ = 4;
2397  *ReplaceSub++ = t[1];
2398  *ReplaceSub++ = t[3];
2399  t += 4;
2400  }
2401 
2402  }
2403  AN.ReplaceScrat[1] = ReplaceSub-AN.ReplaceScrat;
2404  break;
2405 NoRep:
2406  if ( ReplaceType > 0 ) {
2407  C->numrhs = oldtoprhs;
2408  C->Pointer = C->Buffer + oldcpointer;
2409  }
2410  ReplaceType = -1;
2411  t = from;
2412  }
2413 /*
2414  if ( *t == AM.termfunnum && t[1] == FUNHEAD+2
2415  && t[FUNHEAD] == -DOLLAREXPRESSION ) termflag++;
2416 */
2417  if ( *t == DUMMYFUN || *t == DUMMYTEN ) {}
2418  else {
2419  if ( *t < (FUNCTION + WILDOFFSET) ) {
2420  if ( ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2421  || ( functions[*t-FUNCTION].minnumargs > 0 ) )
2422  && ( ( t[2] & DIRTYFLAG ) != 0 ) ) {
2423 /*
2424  Number of arguments is bounded. And we have not checked.
2425 */
2426  WORD *ta = t + FUNHEAD, *tb = t + t[1];
2427  int numarg = 0;
2428  while ( ta < tb ) { numarg++; NEXTARG(ta) }
2429  if ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2430  && ( numarg >= functions[*t-FUNCTION].maxnumargs ) )
2431  goto NormZero;
2432  if ( ( functions[*t-FUNCTION].minnumargs > 0 )
2433  && ( numarg < functions[*t-FUNCTION].minnumargs ) )
2434  goto NormZero;
2435  }
2436 doflags:
2437  if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION].tabl == 0 ) ) {
2438  t[2] &= ~DIRTYFLAG;
2439  t[2] |= DIRTYSYMFLAG;
2440  }
2441  if ( functions[*t-FUNCTION].commute ) { pnco[nnco++] = t; }
2442  else { pcom[ncom++] = t; }
2443  }
2444  else {
2445  if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION-WILDOFFSET].tabl == 0 ) ) {
2446  t[2] &= ~DIRTYFLAG;
2447  t[2] |= DIRTYSYMFLAG;
2448  }
2449  if ( functions[*t-FUNCTION-WILDOFFSET].commute ) {
2450  pnco[nnco++] = t;
2451  }
2452  else { pcom[ncom++] = t; }
2453  }
2454  }
2455 
2456  /* Now hunt for contractible indices */
2457 
2458  if ( ( *t < (FUNCTION + WILDOFFSET)
2459  && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) || (
2460  *t >= (FUNCTION + WILDOFFSET)
2461  && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION ) ) {
2462  if ( *t >= GAMMA && *t <= GAMMASEVEN ) t++;
2463  t += FUNHEAD;
2464  while ( t < r ) {
2465  if ( *t >= AM.OffsetIndex && ( *t >= AM.DumInd
2466  || ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2467  pcon[ncon++] = t;
2468  }
2469  else if ( *t == FUNNYWILD ) { t++; }
2470  t++;
2471  }
2472  }
2473  else {
2474  t += FUNHEAD;
2475  while ( t < r ) {
2476  if ( *t > 0 ) {
2477 /*
2478  Here we should worry about a recursion
2479  A problem is the possibility of a construct
2480  like f(mu+nu)
2481 */
2482  t += *t;
2483  }
2484  else if ( *t <= -FUNCTION ) t++;
2485  else if ( *t == -INDEX ) {
2486  if ( t[1] >= AM.OffsetIndex &&
2487  ( t[1] >= AM.DumInd || ( t[1] < AM.WilInd
2488  && indices[t[1]-AM.OffsetIndex].dimension ) ) )
2489  pcon[ncon++] = t+1;
2490  t += 2;
2491  }
2492  else if ( *t == -SYMBOL ) {
2493  if ( t[1] >= MAXPOWER && t[1] < 2*MAXPOWER ) {
2494  *t = -SNUMBER;
2495  t[1] -= MAXPOWER;
2496  }
2497  else if ( t[1] < -MAXPOWER && t[1] > -2*MAXPOWER ) {
2498  *t = -SNUMBER;
2499  t[1] += MAXPOWER;
2500  }
2501  else t += 2;
2502  }
2503  else t += 2;
2504  }
2505  }
2506  break;
2507  }
2508  t = r;
2509 TryAgain:;
2510  } while ( t < m );
2511  if ( ANsc ) {
2512  AN.cTerm = ANsc;
2513  r = t = ANsr; m = ANsm;
2514  ANsc = ANsm = ANsr = 0;
2515  goto conscan;
2516  }
2517 /*
2518  #] First scan :
2519  #[ Easy denominators :
2520 
2521  Easy denominators are denominators that can be replaced by
2522  negative powers of individual subterms. This may add to all
2523  our sublists.
2524 
2525 */
2526  if ( nden ) {
2527  for ( k = 0, i = 0; i < nden; i++ ) {
2528  t = pden[i];
2529  if ( ( t[2] & DIRTYFLAG ) == 0 ) continue;
2530  r = t + t[1]; m = t + FUNHEAD;
2531  if ( m >= r ) {
2532  for ( j = i+1; j < nden; j++ ) pden[j-1] = pden[j];
2533  nden--;
2534  for ( j = 0; j < nnco; j++ ) if ( pnco[j] == t ) break;
2535  for ( j++; j < nnco; j++ ) pnco[j-1] = pnco[j];
2536  nnco--;
2537  i--;
2538  }
2539  else {
2540  NEXTARG(m);
2541  if ( m >= r ) continue;
2542 /*
2543  We have more than one argument. Split the function.
2544 */
2545  if ( k == 0 ) {
2546  k = 1; to = termout; from = term;
2547  }
2548  while ( from < t ) *to++ = *from++;
2549  m = t + FUNHEAD;
2550  while ( m < r ) {
2551  stop = to;
2552  *to++ = DENOMINATOR;
2553  for ( j = 1; j < FUNHEAD; j++ ) *to++ = 0;
2554  if ( *m < -FUNCTION ) *to++ = *m++;
2555  else if ( *m < 0 ) { *to++ = *m++; *to++ = *m++; }
2556  else {
2557  j = *m; while ( --j >= 0 ) *to++ = *m++;
2558  }
2559  stop[1] = WORDDIF(to,stop);
2560  }
2561  from = r;
2562  if ( i == nden - 1 ) {
2563  stop = term + *term;
2564  while ( from < stop ) *to++ = *from++;
2565  i = *termout = WORDDIF(to,termout);
2566  to = term; from = termout;
2567  while ( --i >= 0 ) *to++ = *from++;
2568  goto Restart;
2569  }
2570  }
2571  }
2572  for ( i = 0; i < nden; i++ ) {
2573  t = pden[i];
2574  if ( ( t[2] & DIRTYFLAG ) == 0 ) continue;
2575  t[2] = 0;
2576  if ( t[FUNHEAD] == -SYMBOL ) {
2577  WORD change;
2578  t += FUNHEAD+1;
2579  change = ExtraSymbol(*t,-1,nsym,ppsym,&ncoef);
2580  nsym += change;
2581  ppsym += change << 1;
2582  goto DropDen;
2583  }
2584  else if ( t[FUNHEAD] == -SNUMBER ) {
2585  t += FUNHEAD+1;
2586  if ( *t == 0 ) goto NormInf;
2587  if ( *t < 0 ) { *AT.WorkPointer = -*t; j = -1; }
2588  else { *AT.WorkPointer = *t; j = 1; }
2589  ncoef = REDLENG(ncoef);
2590  if ( Divvy(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)AT.WorkPointer,j) )
2591  goto FromNorm;
2592  ncoef = INCLENG(ncoef);
2593  goto DropDen;
2594  }
2595  else if ( t[FUNHEAD] == ARGHEAD ) goto NormInf;
2596  else if ( t[FUNHEAD] > 0 && t[FUNHEAD+ARGHEAD] ==
2597  t[FUNHEAD]-ARGHEAD ) {
2598  /* Only one term */
2599  r = t + t[1] - 1;
2600  t += FUNHEAD + ARGHEAD + 1;
2601  j = *r;
2602  m = r - ABS(*r) + 1;
2603  if ( j != 3 || ( ( *m != 1 ) || ( m[1] != 1 ) ) ) {
2604  ncoef = REDLENG(ncoef);
2605  if ( DivRat(BHEAD (UWORD *)AT.n_coef,ncoef,(UWORD *)m,REDLENG(j),(UWORD *)AT.n_coef,&ncoef) ) goto FromNorm;
2606  ncoef = INCLENG(ncoef);
2607  j = ABS(j) - 3;
2608  t[-FUNHEAD-ARGHEAD] -= j;
2609  t[-ARGHEAD-1] -= j;
2610  t[-1] -= j;
2611  m[0] = m[1] = 1;
2612  m[2] = 3;
2613  }
2614  while ( t < m ) {
2615  r = t + t[1];
2616  if ( *t == SYMBOL || *t == DOTPRODUCT ) {
2617  k = t[1];
2618  pden[i][1] -= k;
2619  pden[i][FUNHEAD] -= k;
2620  pden[i][FUNHEAD+ARGHEAD] -= k;
2621  m -= k;
2622  stop = m + 3;
2623  tt = to = t;
2624  from = r;
2625  if ( *t == SYMBOL ) {
2626  t += 2;
2627  while ( t < r ) {
2628  WORD change;
2629  change = ExtraSymbol(*t,-t[1],nsym,ppsym,&ncoef);
2630  nsym += change;
2631  ppsym += change << 1;
2632  t += 2;
2633  }
2634  }
2635  else {
2636  t += 2;
2637  while ( t < r ) {
2638  *ppdot++ = *t++;
2639  *ppdot++ = *t++;
2640  *ppdot++ = -*t++;
2641  ndot++;
2642  }
2643  }
2644  while ( to < stop ) *to++ = *from++;
2645  r = tt;
2646  }
2647  t = r;
2648  }
2649  if ( pden[i][1] == 4+FUNHEAD+ARGHEAD ) {
2650 DropDen:
2651  for ( j = 0; j < nnco; j++ ) {
2652  if ( pden[i] == pnco[j] ) {
2653  --nnco;
2654  while ( j < nnco ) {
2655  pnco[j] = pnco[j+1];
2656  j++;
2657  }
2658  break;
2659  }
2660  }
2661  pden[i--] = pden[--nden];
2662  }
2663  }
2664  }
2665  }
2666 /*
2667  #] Easy denominators :
2668  #[ Index Contractions :
2669 */
2670  if ( ndel ) {
2671  t = pdel;
2672  for ( i = 0; i < ndel; i += 2 ) {
2673  if ( t[0] == t[1] ) {
2674  if ( t[0] == EMPTYINDEX ) {}
2675  else if ( *t < AM.OffsetIndex ) {
2676  k = AC.FixIndices[*t];
2677  if ( k < 0 ) { j = -1; k = -k; }
2678  else if ( k > 0 ) j = 1;
2679  else goto NormZero;
2680  goto WithFix;
2681  }
2682  else if ( *t >= AM.DumInd ) {
2683  k = AC.lDefDim;
2684  if ( k ) goto docontract;
2685  }
2686  else if ( *t >= AM.WilInd ) {
2687  k = indices[*t-AM.OffsetIndex-WILDOFFSET].dimension;
2688  if ( k ) goto docontract;
2689  }
2690  else if ( ( k = indices[*t-AM.OffsetIndex].dimension ) != 0 ) {
2691 docontract:
2692  if ( k > 0 ) {
2693  j = 1;
2694 WithFix: shortnum = k;
2695  ncoef = REDLENG(ncoef);
2696  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)(&shortnum),j) )
2697  goto FromNorm;
2698  ncoef = INCLENG(ncoef);
2699  }
2700  else {
2701  WORD change;
2702  change = ExtraSymbol((WORD)(-k),(WORD)1,nsym,ppsym,&ncoef);
2703  nsym += change;
2704  ppsym += change << 1;
2705  }
2706  t[1] = pdel[ndel-1];
2707  t[0] = pdel[ndel-2];
2708 HaveCon:
2709  ndel -= 2;
2710  i -= 2;
2711  }
2712  }
2713  else {
2714  if ( *t < AM.OffsetIndex && t[1] < AM.OffsetIndex ) goto NormZero;
2715  j = *t - AM.OffsetIndex;
2716  if ( j >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2717  || ( *t < AM.WilInd && indices[j].dimension ) ) ) {
2718  for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2719  if ( *t == *m ) {
2720  *t = m[1];
2721  *m++ = pdel[ndel-2];
2722  *m = pdel[ndel-1];
2723  goto HaveCon;
2724  }
2725  else if ( *t == m[1] ) {
2726  *t = *m;
2727  *m++ = pdel[ndel-2];
2728  *m = pdel[ndel-1];
2729  goto HaveCon;
2730  }
2731  }
2732  }
2733  j = t[1]-AM.OffsetIndex;
2734  if ( j >= 0 && ( ( t[1] >= AM.DumInd && AC.lDefDim )
2735  || ( t[1] < AM.WilInd && indices[j].dimension ) ) ) {
2736  for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2737  if ( t[1] == *m ) {
2738  t[1] = m[1];
2739  *m++ = pdel[ndel-2];
2740  *m = pdel[ndel-1];
2741  goto HaveCon;
2742  }
2743  else if ( t[1] == m[1] ) {
2744  t[1] = *m;
2745  *m++ = pdel[ndel-2];
2746  *m = pdel[ndel-1];
2747  goto HaveCon;
2748  }
2749  }
2750  }
2751  t += 2;
2752  }
2753  }
2754  if ( ndel > 0 ) {
2755  if ( nvec ) {
2756  t = pdel;
2757  for ( i = 0; i < ndel; i++ ) {
2758  if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2759  ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2760  r = pvec + 1;
2761  for ( j = 1; j < nvec; j += 2 ) {
2762  if ( *r == *t ) {
2763  if ( i & 1 ) {
2764  *r = t[-1];
2765  *t-- = pdel[--ndel];
2766  i -= 2;
2767  }
2768  else {
2769  *r = t[1];
2770  t[1] = pdel[--ndel];
2771  i--;
2772  }
2773  *t-- = pdel[--ndel];
2774  break;
2775  }
2776  r += 2;
2777  }
2778  }
2779  t++;
2780  }
2781  }
2782  if ( ndel > 0 && ncon ) {
2783  t = pdel;
2784  for ( i = 0; i < ndel; i++ ) {
2785  if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2786  ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2787  for ( j = 0; j < ncon; j++ ) {
2788  if ( *pcon[j] == *t ) {
2789  if ( i & 1 ) {
2790  *pcon[j] = t[-1];
2791  *t-- = pdel[--ndel];
2792  i -= 2;
2793  }
2794  else {
2795  *pcon[j] = t[1];
2796  t[1] = pdel[--ndel];
2797  i--;
2798  }
2799  *t-- = pdel[--ndel];
2800  didcontr++;
2801  r = pcon[j];
2802  for ( j = 0; j < nnco; j++ ) {
2803  m = pnco[j];
2804  if ( r > m && r < m+m[1] ) {
2805  m[2] |= DIRTYSYMFLAG;
2806  break;
2807  }
2808  }
2809  for ( j = 0; j < ncom; j++ ) {
2810  m = pcom[j];
2811  if ( r > m && r < m+m[1] ) {
2812  m[2] |= DIRTYSYMFLAG;
2813  break;
2814  }
2815  }
2816  for ( j = 0; j < neps; j++ ) {
2817  m = peps[j];
2818  if ( r > m && r < m+m[1] ) {
2819  m[2] |= DIRTYSYMFLAG;
2820  break;
2821  }
2822  }
2823  break;
2824  }
2825  }
2826  }
2827  t++;
2828  }
2829  }
2830  }
2831  }
2832  if ( nvec ) {
2833  t = pvec + 1;
2834  for ( i = 3; i < nvec; i += 2 ) {
2835  k = *t - AM.OffsetIndex;
2836  if ( k >= 0 && ( ( *t > AM.DumInd && AC.lDefDim )
2837  || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2838  r = t + 2;
2839  for ( j = i; j < nvec; j += 2 ) {
2840  if ( *r == *t ) { /* Another dotproduct */
2841  *ppdot++ = t[-1];
2842  *ppdot++ = r[-1];
2843  *ppdot++ = 1;
2844  ndot++;
2845  *r-- = pvec[--nvec];
2846  *r = pvec[--nvec];
2847  *t-- = pvec[--nvec];
2848  *t-- = pvec[--nvec];
2849  i -= 2;
2850  break;
2851  }
2852  r += 2;
2853  }
2854  }
2855  t += 2;
2856  }
2857  if ( nvec > 0 && ncon ) {
2858  t = pvec + 1;
2859  for ( i = 1; i < nvec; i += 2 ) {
2860  k = *t - AM.OffsetIndex;
2861  if ( k >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2862  || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2863  for ( j = 0; j < ncon; j++ ) {
2864  if ( *pcon[j] == *t ) {
2865  *pcon[j] = t[-1];
2866  *t-- = pvec[--nvec];
2867  *t-- = pvec[--nvec];
2868  r = pcon[j];
2869  pcon[j] = pcon[--ncon];
2870  i -= 2;
2871  for ( j = 0; j < nnco; j++ ) {
2872  m = pnco[j];
2873  if ( r > m && r < m+m[1] ) {
2874  m[2] |= DIRTYSYMFLAG;
2875  break;
2876  }
2877  }
2878  for ( j = 0; j < ncom; j++ ) {
2879  m = pcom[j];
2880  if ( r > m && r < m+m[1] ) {
2881  m[2] |= DIRTYSYMFLAG;
2882  break;
2883  }
2884  }
2885  for ( j = 0; j < neps; j++ ) {
2886  m = peps[j];
2887  if ( r > m && r < m+m[1] ) {
2888  m[2] |= DIRTYSYMFLAG;
2889  break;
2890  }
2891  }
2892  break;
2893  }
2894  }
2895  }
2896  t += 2;
2897  }
2898  }
2899  }
2900 /*
2901  #] Index Contractions :
2902  #[ NonCommuting Functions :
2903 */
2904  m = fillsetexp;
2905  if ( nnco ) {
2906  for ( i = 0; i < nnco; i++ ) {
2907  t = pnco[i];
2908  if ( ( *t >= (FUNCTION+WILDOFFSET)
2909  && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
2910  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2911  && functions[*t-FUNCTION].spec == 0 ) ) {
2912  DoRevert(t,m);
2913  if ( didcontr ) {
2914  r = t + FUNHEAD;
2915  t += t[1];
2916  while ( r < t ) {
2917  if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
2918  *r = -SNUMBER;
2919  didcontr--;
2920  pnco[i][2] |= DIRTYSYMFLAG;
2921  }
2922  NEXTARG(r)
2923  }
2924  }
2925  }
2926  }
2927 
2928  /* First should come the code for function properties. */
2929 
2930  /* First we test for symmetric properties and the DIRTYSYMFLAG */
2931 
2932  for ( i = 0; i < nnco; i++ ) {
2933  t = pnco[i];
2934  if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) && *t != DOLLAREXPRESSION ) {
2935  l = 0; /* to make the compiler happy */
2936  if ( ( *t >= (FUNCTION+WILDOFFSET)
2937  && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
2938  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2939  && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
2940  if ( *t >= (FUNCTION+WILDOFFSET) ) {
2941  *t -= WILDOFFSET;
2942  j = FullSymmetrize(BHEAD t,l);
2943  *t += WILDOFFSET;
2944  }
2945  else j = FullSymmetrize(BHEAD t,l);
2946  if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
2947  if ( ( j & 2 ) != 0 ) goto NormZero;
2948  if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
2949  }
2950  }
2951  else t[2] &= ~DIRTYSYMFLAG;
2952  }
2953  }
2954 
2955  /* Non commuting functions are then tested for commutation
2956  rules. If needed their order is exchanged. */
2957 
2958  k = nnco - 1;
2959  for ( i = 0; i < k; i++ ) {
2960  j = i;
2961  while ( Commute(pnco[j],pnco[j+1]) ) {
2962  t = pnco[j]; pnco[j] = pnco[j+1]; pnco[j+1] = t;
2963  l = j-1;
2964  while ( l >= 0 && Commute(pnco[l],pnco[l+1]) ) {
2965  t = pnco[l]; pnco[l] = pnco[l+1]; pnco[l+1] = t;
2966  l--;
2967  }
2968  if ( ++j >= k ) break;
2969  }
2970  }
2971 
2972  /* Finally they are written to output. gamma matrices
2973  are bundled if possible */
2974 
2975  for ( i = 0; i < nnco; i++ ) {
2976  t = pnco[i];
2977  if ( *t >= GAMMA && *t <= GAMMASEVEN ) {
2978  WORD gtype;
2979  to = m;
2980  *m++ = GAMMA;
2981  m++;
2982  FILLFUN(m)
2983  *m++ = stype = t[FUNHEAD]; /* type of string */
2984  j = 0;
2985  nnum = 0;
2986  do {
2987  r = t + t[1];
2988  if ( *t == GAMMAFIVE ) {
2989  gtype = GAMMA5; t += FUNHEAD; goto onegammamatrix; }
2990  else if ( *t == GAMMASIX ) {
2991  gtype = GAMMA6; t += FUNHEAD; goto onegammamatrix; }
2992  else if ( *t == GAMMASEVEN ) {
2993  gtype = GAMMA7; t += FUNHEAD; goto onegammamatrix; }
2994  t += FUNHEAD+1;
2995  while ( t < r ) {
2996  gtype = *t;
2997 onegammamatrix:
2998  if ( gtype == GAMMA5 ) {
2999  if ( j == GAMMA1 ) j = GAMMA5;
3000  else if ( j == GAMMA5 ) j = GAMMA1;
3001  else if ( j == GAMMA7 ) ncoef = -ncoef;
3002  if ( nnum & 1 ) ncoef = -ncoef;
3003  }
3004  else if ( gtype == GAMMA6 || gtype == GAMMA7 ) {
3005  if ( nnum & 1 ) {
3006  if ( gtype == GAMMA6 ) gtype = GAMMA7;
3007  else gtype = GAMMA6;
3008  }
3009  if ( j == GAMMA1 ) j = gtype;
3010  else if ( j == GAMMA5 ) {
3011  j = gtype;
3012  if ( j == GAMMA7 ) ncoef = -ncoef;
3013  }
3014  else if ( j != gtype ) goto NormZero;
3015  else {
3016  shortnum = 2;
3017  ncoef = REDLENG(ncoef);
3018  if ( Mully(BHEAD (UWORD *)AT.n_coef,&ncoef,(UWORD *)(&shortnum),1) ) goto FromNorm;
3019  ncoef = INCLENG(ncoef);
3020  }
3021  }
3022  else {
3023  *m++ = gtype; nnum++;
3024  }
3025  t++;
3026  }
3027 
3028  } while ( ( ++i < nnco ) && ( *(t = pnco[i]) >= GAMMA
3029  && *t <= GAMMASEVEN ) && ( t[FUNHEAD] == stype ) );
3030  i--;
3031  if ( j ) {
3032  k = WORDDIF(m,to) - FUNHEAD-1;
3033  r = m;
3034  from = m++;
3035  while ( --k >= 0 ) *from-- = *--r;
3036  *from = j;
3037  }
3038  to[1] = WORDDIF(m,to);
3039  }
3040  else if ( *t < 0 ) {
3041  *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3042  FILLFUN3(m)
3043  }
3044  else {
3045  k = t[1];
3046  NCOPY(m,t,k);
3047  }
3048  }
3049 
3050  }
3051 /*
3052  #] NonCommuting Functions :
3053  #[ Commuting Functions :
3054 */
3055  if ( ncom ) {
3056  for ( i = 0; i < ncom; i++ ) {
3057  t = pcom[i];
3058  if ( ( *t >= (FUNCTION+WILDOFFSET)
3059  && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
3060  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3061  && functions[*t-FUNCTION].spec == 0 ) ) {
3062  DoRevert(t,m);
3063  if ( didcontr ) {
3064  r = t + FUNHEAD;
3065  t += t[1];
3066  while ( r < t ) {
3067  if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3068  *r = -SNUMBER;
3069  didcontr--;
3070  pcom[i][2] |= DIRTYSYMFLAG;
3071  }
3072  NEXTARG(r)
3073  }
3074  }
3075  }
3076  }
3077 
3078  /* Now we test for symmetric properties and the DIRTYSYMFLAG */
3079 
3080  for ( i = 0; i < ncom; i++ ) {
3081  t = pcom[i];
3082  if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) ) {
3083  l = 0; /* to make the compiler happy */
3084  if ( ( *t >= (FUNCTION+WILDOFFSET)
3085  && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3086  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3087  && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3088  if ( *t >= (FUNCTION+WILDOFFSET) ) {
3089  *t -= WILDOFFSET;
3090  j = FullSymmetrize(BHEAD t,l);
3091  *t += WILDOFFSET;
3092  }
3093  else j = FullSymmetrize(BHEAD t,l);
3094  if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3095  if ( ( j & 2 ) != 0 ) goto NormZero;
3096  if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3097  }
3098  }
3099  else t[2] &= ~DIRTYSYMFLAG;
3100  }
3101  }
3102 /*
3103  Sort the functions
3104  From a purists point of view this can be improved.
3105  There arel slow and fast arguments and no conversions are
3106  taken into account here.
3107 */
3108  for ( i = 1; i < ncom; i++ ) {
3109  for ( j = i; j > 0; j-- ) {
3110  WORD jj,kk;
3111  jj = j-1;
3112  t = pcom[jj];
3113  r = pcom[j];
3114  if ( *t < 0 ) {
3115  if ( *r < 0 ) { if ( *t >= *r ) goto NextI; }
3116  else { if ( -*t <= *r ) goto NextI; }
3117  goto jexch;
3118  }
3119  else if ( *r < 0 ) {
3120  if ( *t < -*r ) goto NextI;
3121  goto jexch;
3122  }
3123  else if ( *t != *r ) {
3124  if ( *t < *r ) goto NextI;
3125 jexch: t = pcom[j]; pcom[j] = pcom[jj]; pcom[jj] = t;
3126  continue;
3127  }
3128  if ( AC.properorderflag ) {
3129  if ( ( *t >= (FUNCTION+WILDOFFSET)
3130  && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
3131  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3132  && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) ) {}
3133  else {
3134  WORD *s1, *s2, *ss1, *ss2;
3135  s1 = t+FUNHEAD; s2 = r+FUNHEAD;
3136  ss1 = t + t[1]; ss2 = r + r[1];
3137  while ( s1 < ss1 && s2 < ss2 ) {
3138  k = CompArg(s1,s2);
3139  if ( k > 0 ) goto jexch;
3140  if ( k < 0 ) goto NextI;
3141  NEXTARG(s1)
3142  NEXTARG(s2)
3143  }
3144  if ( s1 < ss1 ) goto jexch;
3145  goto NextI;
3146  }
3147  k = t[1] - FUNHEAD;
3148  kk = r[1] - FUNHEAD;
3149  t += FUNHEAD;
3150  r += FUNHEAD;
3151  while ( k > 0 && kk > 0 ) {
3152  if ( *t < *r ) goto NextI;
3153  else if ( *t++ > *r++ ) goto jexch;
3154  k--; kk--;
3155  }
3156  if ( k > 0 ) goto jexch;
3157  goto NextI;
3158  }
3159  else
3160  {
3161  k = t[1] - FUNHEAD;
3162  kk = r[1] - FUNHEAD;
3163  t += FUNHEAD;
3164  r += FUNHEAD;
3165  while ( k > 0 && kk > 0 ) {
3166  if ( *t < *r ) goto NextI;
3167  else if ( *t++ > *r++ ) goto jexch;
3168  k--; kk--;
3169  }
3170  if ( k > 0 ) goto jexch;
3171  goto NextI;
3172  }
3173  }
3174 NextI:;
3175  }
3176  for ( i = 0; i < ncom; i++ ) {
3177  t = pcom[i];
3178  if ( *t == THETA || *t == THETA2 ) {
3179  if ( ( k = DoTheta(BHEAD t) ) == 0 ) goto NormZero;
3180  else if ( k < 0 ) {
3181  k = t[1];
3182  NCOPY(m,t,k);
3183  }
3184  }
3185  else if ( *t == DELTA2 || *t == DELTAP ) {
3186  if ( ( k = DoDelta(t) ) == 0 ) goto NormZero;
3187  else if ( k < 0 ) {
3188  k = t[1];
3189  NCOPY(m,t,k);
3190  }
3191  }
3192  else if ( *t == AR.PolyFun ) {
3193  if ( AR.PolyFunType == 1 ) { /* Regular PolyFun with one argument */
3194  if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3195  t[1] == FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) goto NormZero;
3196  if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3197  k = t[1];
3198  NCOPY(m,t,k);
3199  }
3200  else if ( AR.PolyFunType == 2 ) { /* PolyRatFun. Two arguments */
3201 /*
3202  First check for zeroes.
3203 */
3204  if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3205  t[1] > FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) {
3206  u = t + FUNHEAD + 2;
3207  if ( *u < 0 ) {
3208  if ( *u <= -FUNCTION ) {}
3209  else if ( t[1] == FUNHEAD+4 && t[FUNHEAD+2] == -SNUMBER
3210  && t[FUNHEAD+3] == 0 ) goto NormPRF;
3211  else if ( t[1] == FUNHEAD+4 ) goto NormZero;
3212  }
3213  else if ( t[1] == *u+FUNHEAD+2 ) goto NormZero;
3214  }
3215  if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3216  k = t[1];
3217  NCOPY(m,t,k);
3218  }
3219  }
3220  else if ( *t > 0 ) {
3221  k = t[1];
3222  NCOPY(m,t,k);
3223  }
3224  else {
3225  *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3226  FILLFUN3(m)
3227  }
3228  }
3229  }
3230 /*
3231  #] Commuting Functions :
3232  #[ LeviCivita tensors :
3233 */
3234  if ( neps ) {
3235  to = m;
3236  for ( i = 0; i < neps; i++ ) { /* Put the indices in order */
3237  t = peps[i];
3238  if ( ( t[2] & DIRTYSYMFLAG ) != DIRTYSYMFLAG ) continue;
3239  t[2] &= ~DIRTYSYMFLAG;
3240  if ( AR.Eside == LHSIDE || AR.Eside == LHSIDEX ) {
3241  /* Potential problems with FUNNYWILD */
3242 /*
3243  First make sure all FUNNIES are at the end.
3244  Then sort separately
3245 */
3246  r = t + FUNHEAD;
3247  m = tt = t + t[1];
3248  while ( r < m ) {
3249  if ( *r != FUNNYWILD ) { r++; continue; }
3250  k = r[1]; u = r + 2;
3251  while ( u < tt ) {
3252  u[-2] = *u;
3253  if ( *u != FUNNYWILD ) ncoef = -ncoef;
3254  u++;
3255  }
3256  tt[-2] = FUNNYWILD; tt[-1] = k; m -= 2;
3257  }
3258  t += FUNHEAD;
3259  do {
3260  for ( r = t + 1; r < m; r++ ) {
3261  if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3262  else if ( *r == *t ) goto NormZero;
3263  }
3264  t++;
3265  } while ( t < m );
3266  do {
3267  for ( r = t + 2; r < tt; r += 2 ) {
3268  if ( r[1] < t[1] ) {
3269  k = r[1]; r[1] = t[1]; t[1] = k; ncoef = -ncoef; }
3270  else if ( r[1] == t[1] ) goto NormZero;
3271  }
3272  t += 2;
3273  } while ( t < tt );
3274  }
3275  else {
3276  m = t + t[1];
3277  t += FUNHEAD;
3278  do {
3279  for ( r = t + 1; r < m; r++ ) {
3280  if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3281  else if ( *r == *t ) goto NormZero;
3282  }
3283  t++;
3284  } while ( t < m );
3285  }
3286  }
3287 
3288  /* Sort the tensors */
3289 
3290  for ( i = 0; i < (neps-1); i++ ) {
3291  t = peps[i];
3292  for ( j = i+1; j < neps; j++ ) {
3293  r = peps[j];
3294  if ( t[1] > r[1] ) {
3295  peps[i] = m = r; peps[j] = r = t; t = m;
3296  }
3297  else if ( t[1] == r[1] ) {
3298  k = t[1] - FUNHEAD;
3299  m = t + FUNHEAD;
3300  r += FUNHEAD;
3301  do {
3302  if ( *r < *m ) {
3303  m = peps[j]; peps[j] = t; peps[i] = t = m;
3304  break;
3305  }
3306  else if ( *r++ > *m++ ) break;
3307  } while ( --k > 0 );
3308  }
3309  }
3310  }
3311  m = to;
3312  for ( i = 0; i < neps; i++ ) {
3313  t = peps[i];
3314  k = t[1];
3315  NCOPY(m,t,k);
3316  }
3317  }
3318 /*
3319  #] LeviCivita tensors :
3320  #[ Delta :
3321 */
3322  if ( ndel ) {
3323  r = t = pdel;
3324  for ( i = 0; i < ndel; i += 2, r += 2 ) {
3325  if ( r[1] < r[0] ) { k = *r; *r = r[1]; r[1] = k; }
3326  }
3327  for ( i = 2; i < ndel; i += 2, t += 2 ) {
3328  r = t + 2;
3329  for ( j = i; j < ndel; j += 2 ) {
3330  if ( *r > *t ) { r += 2; }
3331  else if ( *r < *t ) {
3332  k = *r; *r++ = *t; *t++ = k;
3333  k = *r; *r++ = *t; *t-- = k;
3334  }
3335  else {
3336  if ( *++r < t[1] ) {
3337  k = *r; *r = t[1]; t[1] = k;
3338  }
3339  r++;
3340  }
3341  }
3342  }
3343  t = pdel;
3344  *m++ = DELTA;
3345  *m++ = ndel + 2;
3346  i = ndel;
3347  NCOPY(m,t,i);
3348  }
3349 /*
3350  #] Delta :
3351  #[ Loose Vectors/Indices :
3352 */
3353  if ( nind ) {
3354  t = pind;
3355  for ( i = 0; i < nind; i++ ) {
3356  r = t + 1;
3357  for ( j = i+1; j < nind; j++ ) {
3358  if ( *r < *t ) {
3359  k = *r; *r = *t; *t = k;
3360  }
3361  r++;
3362  }
3363  t++;
3364  }
3365  t = pind;
3366  *m++ = INDEX;
3367  *m++ = nind + 2;
3368  i = nind;
3369  NCOPY(m,t,i);
3370  }
3371 /*
3372  #] Loose Vectors/Indices :
3373  #[ Vectors :
3374 */
3375  if ( nvec ) {
3376  t = pvec;
3377  for ( i = 2; i < nvec; i += 2 ) {
3378  r = t + 2;
3379  for ( j = i; j < nvec; j += 2 ) {
3380  if ( *r == *t ) {
3381  if ( *++r < t[1] ) {
3382  k = *r; *r = t[1]; t[1] = k;
3383  }
3384  r++;
3385  }
3386  else if ( *r < *t ) {
3387  k = *r; *r++ = *t; *t++ = k;
3388  k = *r; *r++ = *t; *t-- = k;
3389  }
3390  else { r += 2; }
3391  }
3392  t += 2;
3393  }
3394  t = pvec;
3395  *m++ = VECTOR;
3396  *m++ = nvec + 2;
3397  i = nvec;
3398  NCOPY(m,t,i);
3399  }
3400 /*
3401  #] Vectors :
3402  #[ Dotproducts :
3403 */
3404  if ( ndot ) {
3405  to = m;
3406  m = t = pdot;
3407  i = ndot;
3408  while ( --i >= 0 ) {
3409  if ( *t > t[1] ) { j = *t; *t = t[1]; t[1] = j; }
3410  t += 3;
3411  }
3412  t = m;
3413  ndot *= 3;
3414  m += ndot;
3415  while ( t < (m-3) ) {
3416  r = t + 3;
3417  do {
3418  if ( *r == *t ) {
3419  if ( *++r == *++t ) {
3420  r++;
3421  if ( ( *r < MAXPOWER && t[1] < MAXPOWER )
3422  || ( *r > -MAXPOWER && t[1] > -MAXPOWER ) ) {
3423  t++;
3424  *t += *r;
3425  if ( *t > MAXPOWER || *t < -MAXPOWER ) {
3426  MLOCK(ErrorMessageLock);
3427  MesPrint("Exponent of dotproduct out of range: %d",*t);
3428  MUNLOCK(ErrorMessageLock);
3429  goto NormMin;
3430  }
3431  ndot -= 3;
3432  *r-- = *--m;
3433  *r-- = *--m;
3434  *r = *--m;
3435  if ( !*t ) {
3436  ndot -= 3;
3437  *t-- = *--m;
3438  *t-- = *--m;
3439  *t = *--m;
3440  t -= 3;
3441  break;
3442  }
3443  }
3444  else if ( *r < *++t ) {
3445  k = *r; *r++ = *t; *t = k;
3446  }
3447  else r++;
3448  t -= 2;
3449  }
3450  else if ( *r < *t ) {
3451  k = *r; *r++ = *t; *t++ = k;
3452  k = *r; *r++ = *t; *t = k;
3453  t -= 2;
3454  }
3455  else { r += 2; t--; }
3456  }
3457  else if ( *r < *t ) {
3458  k = *r; *r++ = *t; *t++ = k;
3459  k = *r; *r++ = *t; *t++ = k;
3460  k = *r; *r++ = *t; *t = k;
3461  t -= 2;
3462  }
3463  else { r += 3; }
3464  } while ( r < m );
3465  t += 3;
3466  }
3467  m = to;
3468  t = pdot;
3469  if ( ( i = ndot ) > 0 ) {
3470  *m++ = DOTPRODUCT;
3471  *m++ = i + 2;
3472  NCOPY(m,t,i);
3473  }
3474  }
3475 /*
3476  #] Dotproducts :
3477  #[ Symbols :
3478 */
3479  if ( nsym ) {
3480  nsym <<= 1;
3481  t = psym;
3482  *m++ = SYMBOL;
3483  r = m;
3484  *m++ = ( i = nsym ) + 2;
3485  if ( i ) { do {
3486  if ( !*t ) {
3487  if ( t[1] < (2*MAXPOWER) ) { /* powers of i */
3488  if ( t[1] & 1 ) { *m++ = 0; *m++ = 1; }
3489  else *r -= 2;
3490  if ( *++t & 2 ) ncoef = -ncoef;
3491  t++;
3492  }
3493  }
3494  else if ( *t <= NumSymbols && *t > -2*MAXPOWER ) { /* Put powers in range */
3495  if ( ( ( ( t[1] > symbols[*t].maxpower ) && ( symbols[*t].maxpower < MAXPOWER ) ) ||
3496  ( ( t[1] < symbols[*t].minpower ) && ( symbols[*t].minpower > -MAXPOWER ) ) ) &&
3497  ( t[1] < 2*MAXPOWER ) && ( t[1] > -2*MAXPOWER ) ) {
3498  if ( i <= 2 || t[2] != *t ) goto NormZero;
3499  }
3500  if ( AN.ncmod == 1 && ( AC.modmode & ALSOPOWERS ) != 0 ) {
3501  if ( AC.cmod[0] == 1 ) t[1] = 0;
3502  else if ( t[1] >= 0 ) t[1] = 1 + (t[1]-1)%(AC.cmod[0]-1);
3503  else {
3504  t[1] = -1 - (-t[1]-1)%(AC.cmod[0]-1);
3505  if ( t[1] < 0 ) t[1] += (AC.cmod[0]-1);
3506  }
3507  }
3508  if ( ( t[1] < (2*MAXPOWER) && t[1] >= MAXPOWER )
3509  || ( t[1] > -(2*MAXPOWER) && t[1] <= -MAXPOWER ) ) {
3510  MLOCK(ErrorMessageLock);
3511  MesPrint("Exponent out of range: %d",t[1]);
3512  MUNLOCK(ErrorMessageLock);
3513  goto NormMin;
3514  }
3515  if ( t[1] ) {
3516  *m++ = *t++;
3517  *m++ = *t++;
3518  }
3519  else { *r -= 2; t += 2; }
3520  }
3521  else {
3522  *m++ = *t++; *m++ = *t++;
3523  }
3524  } while ( (i-=2) > 0 ); }
3525  if ( *r <= 2 ) m = r-1;
3526  }
3527 /*
3528  #] Symbols :
3529  #[ Errors and Finish :
3530 */
3531  stop = (WORD *)(((UBYTE *)(termout)) + AM.MaxTer);
3532  i = ABS(ncoef);
3533  if ( ( m + i ) > stop ) {
3534  MLOCK(ErrorMessageLock);
3535  MesPrint("Term too complex during normalization");
3536  MUNLOCK(ErrorMessageLock);
3537  goto NormMin;
3538  }
3539  if ( ReplaceType >= 0 ) {
3540  t = AT.n_coef;
3541  i--;
3542  NCOPY(m,t,i);
3543  *m++ = ncoef;
3544  t = termout;
3545  *t = WORDDIF(m,t);
3546  if ( ReplaceType == 0 ) {
3547  AT.WorkPointer = termout+*termout;
3548  WildFill(BHEAD term,termout,AN.ReplaceScrat);
3549  termout = term + *term;
3550  }
3551  else {
3552  AT.WorkPointer = r = termout + *termout;
3553  WildFill(BHEAD r,termout,AN.ReplaceScrat);
3554  i = *r; m = term;
3555  NCOPY(m,r,i);
3556  termout = m;
3557 
3558 
3559  r = m = term;
3560  r += *term; r -= ABS(r[-1]);
3561  m++;
3562  while ( m < r ) {
3563  if ( *m > FUNCTION && m[1] > FUNHEAD &&
3564  functions[*m-FUNCTION].spec != TENSORFUNCTION )
3565  m[2] |= DIRTYFLAG;
3566  m += m[1];
3567  }
3568  }
3569 /*
3570  #[ normalize replacements :
3571 
3572  At this point we have the problem that we may have to treat functions
3573  with a dirtyflag. In the original setting DIRTYFLAG is replaced
3574  in the redo by DIRTYSYMFLAG but that doesn't take care of things
3575  like f(y*x) -> f(x*y) etc. This is why we need to redo the arguments
3576  in the style that TestSub uses for dirty arguments, But this again
3577  involves calls to Normalize itself (and to the sorting system).
3578  This is the reason why Normalize has to be reentrant.
3579 
3580  To do things 100% right we have to call TestSub and potentially
3581  invoke Generator, because there could be Table elements that
3582  suddenly can be substituted!
3583 
3584  It seems best to get the term through Generator again, but that means
3585  we have to catch the output via the sort mechanism. It may be a bit
3586  wasteful, but it is definitely the most correct.
3587 */
3588 #ifdef OLDNORMREPLACE
3589  {
3590  WORD oldpolynorm = AN.PolyNormFlag;
3591  WORD *oldcterm = AN.cTerm, *tstop, *argstop, *rnext, *tt, *wt;
3592  WORD *oldwork, olddefer;
3593  LONG newspace, oldspace, numterms;
3594  AT.WorkPointer = termout;
3595  if ( AT.WorkPointer < term + *term && term >= AT.WorkSpace
3596  && term < AT.WorkTop ) AT.WorkPointer = term + *term;
3597 /*
3598  To do things 100% right we have to call TestSub and potentially
3599  invoke Generator, because there could be Table elements that
3600  suddenly can be substituted!
3601  That last possibility we will omit!
3602 */
3603  tstop = term + *term; tstop -= ABS(tstop[-1]);
3604  t = term +1;
3605  while ( t < tstop ) {
3606  if ( *t >= FUNCTION && ( ( t[2] & DIRTYFLAG ) != 0 )
3607  && ( functions[*t-FUNCTION].spec == 0 ) ) {
3608  VOID *oldcompareroutine = AR.CompareRoutine;
3609  r = t + FUNHEAD; argstop = t + t[1];
3610  while ( r < argstop ) {
3611  if ( *r > 0 && ( r[1] != 0 ) ) {
3612  m = r + ARGHEAD; rnext = r + *r;
3613  oldwork = AT.WorkPointer;
3614  olddefer = AR.DeferFlag;
3615  AR.DeferFlag = 0;
3616  if ( *t == AR.PolyFun && AR.PolyFunType == 2 ) {
3617  AR.CompareRoutine = &CompareSymbols;
3618  }
3619  NewSort(BHEAD0);
3620  m = r + ARGHEAD; rnext = r + *r;
3621  while ( m < rnext ) {
3622  i = *m; tt = m; wt = oldwork;
3623  NCOPY(wt,tt,i);
3624  AT.WorkPointer = wt;
3625  if ( Generator(BHEAD oldwork,AR.Cnumlhs) ) {
3626  LowerSortLevel(); goto FromNorm;
3627  }
3628  m += *m;
3629  }
3630  AT.WorkPointer = (WORD *)(((UBYTE *)(oldwork)) + AM.MaxTer);
3631  if ( AT.WorkPointer > AT.WorkTop ) goto OverWork;
3632  m = AT.WorkPointer;
3633  if ( EndSort(BHEAD m,1) < 0 ) goto FromNorm;
3634  if ( *t == AR.PolyFun && AR.PolyFunType == 2 ) {
3635  AR.CompareRoutine = oldcompareroutine;
3636  }
3637 /*
3638  Now we have to analyse the output
3639  Count terms and space
3640 */
3641  AR.DeferFlag = olddefer;
3642  numterms = 0; wt = m;
3643  while ( *wt ) { numterms++; wt += *wt; }
3644  newspace = wt - m;
3645  oldspace = *r - ARGHEAD;
3646 /*
3647  Now the special cases
3648 */
3649  if ( numterms == 0 ) {
3650  m[0] = -SNUMBER; m[1] = 0;
3651  newspace = 2;
3652  }
3653  else if ( numterms == 1 ) {
3654  if ( *m == 4+FUNHEAD && m[3+FUNHEAD] == 3
3655  && m[2+FUNHEAD] == 1 && m[1+FUNHEAD] == 1
3656  && m[1] >= FUNCTION ) {
3657  m[0] = -m[1];
3658  newspace = 1;
3659  }
3660  else if ( *m == 8 && m[7] == 3
3661  && m[6] == 1 && m[5] == 1
3662  && m[1] == SYMBOL && m[4] == 1 ) {
3663  m[0] = -SYMBOL; m[1] = m[3];
3664  newspace = 2;
3665  }
3666  else if ( *m == 7 && m[6] == 3
3667  && m[5] == 1 && m[4] == 1
3668  && m[1] == INDEX ) {
3669  if ( m[3] >= 0 ) {
3670  m[0] = -INDEX; m[1] = m[3];
3671  }
3672  else {
3673  m[0] = -VECTOR; m[1] = m[3];
3674  }
3675  newspace = 2;
3676  }
3677  else if ( *m == 7 && m[6] == -3
3678  && m[5] == 1 && m[4] == 1
3679  && m[1] == INDEX && m[3] < 0 ) {
3680  m[0] = -MINVECTOR; m[1] = m[3];
3681  newspace = 2;
3682  }
3683  else if ( *m == 4
3684  && m[2] == 1 && (UWORD)(m[1]) <= MAXPOSITIVE ) {
3685  m[0] = -SNUMBER;
3686  if ( m[3] < 0 ) m[1] = -m[1];
3687  newspace = 2;
3688  }
3689  }
3690 /*
3691  Now the old argument takes oldspace spaces.
3692  The new argument takes newspace places.
3693  The new argument sits at m. There should be enough
3694  space in the term to accept it, but we may have to
3695  move the tail of the term
3696 */
3697  if ( newspace <= 2 ) {
3698  oldspace = *r;
3699  i = oldspace - newspace;
3700  *r = *m;
3701  if ( newspace > 1 ) r[1] = m[1];
3702  m = r + newspace;
3703  wt = rnext;
3704  tt = term + *term;
3705  while ( wt < tt ) *m++ = *wt++;
3706  *term -= i;
3707  t[1] -= i;
3708  tstop -= i;
3709  argstop -= i;
3710  }
3711  else if ( oldspace == newspace ) {
3712  i = newspace; tt = r+ARGHEAD; wt = m;
3713  NCOPY(tt,wt,i);
3714  r[1] = 0;
3715  }
3716  else if ( oldspace > newspace ) {
3717  i = newspace; tt = r+ARGHEAD; wt = m;
3718  NCOPY(tt,wt,i);
3719  wt = rnext; m = term + *term;
3720  while ( wt < m ) *tt++ = *wt++;
3721  i = oldspace - newspace;
3722  *term -= i;
3723  t[1] -= i;
3724  tstop -= i;
3725  argstop -= i;
3726  *r -= i;
3727  r[1] = 0;
3728  }
3729  else if ( (*term+newspace-oldspace)*sizeof(WORD) > AM.MaxTer ) {
3730  MLOCK(ErrorMessageLock);
3731  MesPrint("Term too complex. Maybe increasing MaxTermSize can help");
3732  MesCall("Norm");
3733  MUNLOCK(ErrorMessageLock);
3734  return(-1);
3735  }
3736  else {
3737  i = newspace - oldspace;
3738  tt = term + *term; wt = rnext;
3739  while ( tt > rnext ) { tt--; tt[i] = tt[0]; }
3740  *term += i;
3741  t[1] += i;
3742  tstop += i;
3743  argstop += i;
3744  *r += i;
3745  i = newspace; tt = r+ARGHEAD; wt = m;
3746  NCOPY(tt,wt,i);
3747  r[1] = 0;
3748  }
3749  AT.WorkPointer = oldwork;
3750  }
3751  NEXTARG(r)
3752  }
3753  }
3754  if ( *t >= FUNCTION && ( t[2] & DIRTYFLAG ) != 0 ) {
3755  t[2] &= ~DIRTYFLAG;
3756  if ( functions[*t-FUNCTION].symmetric ) t[2] |= DIRTYSYMFLAG;
3757  }
3758  t += t[1];
3759  }
3760 
3761  AN.PolyNormFlag = oldpolynorm;
3762  AN.cTerm = oldcterm;
3763  }
3764  {
3765  WORD *oldwork = AT.WorkPointer;
3766  WORD olddefer = AR.DeferFlag;
3767  AR.DeferFlag = 0;
3768  NewSort(BHEAD0);
3769  if ( Generator(BHEAD term,AR.Cnumlhs) ) {
3770  LowerSortLevel(); goto FromNorm;
3771  }
3772  AT.WorkPointer = oldwork;
3773  if ( EndSort(BHEAD term,1) < 0 ) goto FromNorm;
3774  if ( *term == 0 ) goto NormZero;
3775  AR.DeferFlag = olddefer;
3776  }
3777 #endif
3778 /*
3779  #] normalize replacements :
3780 */
3781 #ifdef OLDNORMREPLACE
3782  AT.WorkPointer = termout;
3783  if ( ReplaceType == 0 ) {
3784  regval = 1;
3785  goto Restart;
3786  }
3787 #endif
3788 /*
3789  The next 'reset' cannot be done. We still need the expression
3790  in the buffer. Note though that this may cause a runaway pointer
3791  if we are not very careful.
3792 
3793  C->numrhs = oldtoprhs;
3794  C->Pointer = C->Buffer + oldcpointer;
3795 */
3796  return(1);
3797  }
3798  else {
3799  t = termout;
3800  k = WORDDIF(m,t);
3801  *t = k + i;
3802  m = term;
3803  NCOPY(m,t,k);
3804  i--;
3805  t = AT.n_coef;
3806  NCOPY(m,t,i);
3807  *m++ = ncoef;
3808  }
3809 RegEnd:
3810  AT.WorkPointer = termout;
3811  if ( termout < term + *term && termout >= term ) AT.WorkPointer = term + *term;
3812 /*
3813  if ( termflag ) { We have to assign the term to $variable(s)
3814  TermAssign(term);
3815  }
3816 */
3817  return(regval);
3818 
3819 NormInf:
3820  MLOCK(ErrorMessageLock);
3821  MesPrint("Division by zero during normalization");
3822  MUNLOCK(ErrorMessageLock);
3823  goto NormMin;
3824 
3825 NormZZ:
3826  MLOCK(ErrorMessageLock);
3827  MesPrint("0^0 during normalization of term");
3828  MUNLOCK(ErrorMessageLock);
3829  goto NormMin;
3830 
3831 NormPRF:
3832  MLOCK(ErrorMessageLock);
3833  MesPrint("0/0 in polyratfun during normalization of term");
3834  MUNLOCK(ErrorMessageLock);
3835  goto NormMin;
3836 
3837 NormZero:
3838  *term = 0;
3839  AT.WorkPointer = termout;
3840  return(regval);
3841 
3842 NormMin:
3843  return(-1);
3844 
3845 FromNorm:
3846  MLOCK(ErrorMessageLock);
3847  MesCall("Norm");
3848  MUNLOCK(ErrorMessageLock);
3849  return(-1);
3850 
3851 #ifdef OLDNORMREPLACE
3852 OverWork:
3853  MLOCK(ErrorMessageLock);
3854  MesWork();
3855  MesCall("Norm");
3856  MUNLOCK(ErrorMessageLock);
3857  return(-1);
3858 #endif
3859 
3860 /*
3861  #] Errors and Finish :
3862 */
3863 }
3864 
3865 /*
3866  #] Normalize :
3867  #[ ExtraSymbol :
3868 */
3869 
3870 WORD ExtraSymbol(WORD sym, WORD pow, WORD nsym, WORD *ppsym, WORD *ncoef)
3871 {
3872  WORD *m, i;
3873  i = nsym;
3874  m = ppsym - 2;
3875  while ( i > 0 ) {
3876  if ( sym == *m ) {
3877  m++;
3878  if ( pow > 2*MAXPOWER || pow < -2*MAXPOWER
3879  || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
3880  MLOCK(ErrorMessageLock);
3881  MesPrint("Illegal wildcard power combination.");
3882  MUNLOCK(ErrorMessageLock);
3883  Terminate(-1);
3884  }
3885  *m += pow;
3886 
3887  if ( ( sym <= NumSymbols && sym > -MAXPOWER )
3888  && ( symbols[sym].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
3889  *m %= symbols[sym].maxpower;
3890  if ( *m < 0 ) *m += symbols[sym].maxpower;
3891  if ( ( symbols[sym].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
3892  if ( ( ( symbols[sym].maxpower & 1 ) == 0 ) &&
3893  ( *m >= symbols[sym].maxpower/2 ) ) {
3894  *m -= symbols[sym].maxpower/2; *ncoef = -*ncoef;
3895  }
3896  }
3897  }
3898 
3899  if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
3900  MLOCK(ErrorMessageLock);
3901  MesPrint("Power overflow during normalization");
3902  MUNLOCK(ErrorMessageLock);
3903  return(-1);
3904  }
3905  if ( !*m ) {
3906  m--;
3907  while ( i < nsym )
3908  { *m = m[2]; m++; *m = m[2]; m++; i++; }
3909  return(-1);
3910  }
3911  return(0);
3912  }
3913  else if ( sym < *m ) {
3914  m -= 2;
3915  i--;
3916  }
3917  else break;
3918  }
3919  m = ppsym;
3920  while ( i < nsym )
3921  { m--; m[2] = *m; m--; m[2] = *m; i++; }
3922  *m++ = sym;
3923  *m = pow;
3924  return(1);
3925 }
3926 
3927 /*
3928  #] ExtraSymbol :
3929  #[ DoTheta :
3930 */
3931 
3932 WORD DoTheta(PHEAD WORD *t)
3933 {
3934  GETBIDENTITY
3935  WORD k, *r1, *r2, *tstop, type;
3936  WORD ia, *ta, *tb, *stopa, *stopb;
3937  if ( AC.BracketNormalize ) return(-1);
3938  type = *t;
3939  k = t[1];
3940  tstop = t + k;
3941  t += FUNHEAD;
3942  if ( k <= FUNHEAD ) return(1);
3943  r1 = t;
3944  NEXTARG(r1)
3945  if ( r1 == tstop ) {
3946 /*
3947  One argument
3948 */
3949  if ( *t == ARGHEAD ) {
3950  if ( type == THETA ) return(1);
3951  else return(0); /* THETA2 */
3952  }
3953  if ( *t < 0 ) {
3954  if ( *t == -SNUMBER ) {
3955  if ( t[1] < 0 ) return(0);
3956  else {
3957  if ( type == THETA2 && t[1] == 0 ) return(0);
3958  else return(1);
3959  }
3960  }
3961  return(-1);
3962  }
3963  k = t[*t-1];
3964  if ( *t == ABS(k)+1+ARGHEAD ) {
3965  if ( k > 0 ) return(1);
3966  else return(0);
3967  }
3968  return(-1);
3969  }
3970 /*
3971  At least two arguments
3972 */
3973  r2 = r1;
3974  NEXTARG(r2)
3975  if ( r2 < tstop ) return(-1); /* More than 2 arguments */
3976 /*
3977  Note now that zero has to be treated specially
3978  We take the criteria from the symmetrize routine
3979 */
3980  if ( *t == -SNUMBER && *r1 == -SNUMBER ) {
3981  if ( t[1] > r1[1] ) return(0);
3982  else if ( t[1] < r1[1] ) {
3983  return(1);
3984  }
3985  else if ( type == THETA ) return(1);
3986  else return(0); /* THETA2 */
3987  }
3988  else if ( t[1] == 0 && *t == -SNUMBER ) {
3989  if ( *r1 > 0 ) { }
3990  else if ( *t < *r1 ) return(1);
3991  else if ( *t > *r1 ) return(0);
3992  }
3993  else if ( r1[1] == 0 && *r1 == -SNUMBER ) {
3994  if ( *t > 0 ) { }
3995  else if ( *t < *r1 ) return(1);
3996  else if ( *t > *r1 ) return(0);
3997  }
3998  r2 = AT.WorkPointer;
3999  if ( *t < 0 ) {
4000  ta = r2;
4001  ToGeneral(t,ta,0);
4002  r2 += *r2;
4003  }
4004  else ta = t;
4005  if ( *r1 < 0 ) {
4006  tb = r2;
4007  ToGeneral(r1,tb,0);
4008  }
4009  else tb = r1;
4010  stopa = ta + *ta;
4011  stopb = tb + *tb;
4012  ta += ARGHEAD; tb += ARGHEAD;
4013  while ( ta < stopa ) {
4014  if ( tb >= stopb ) return(0);
4015  if ( ( ia = CompareTerms(BHEAD ta,tb,(WORD)1) ) < 0 ) return(0);
4016  if ( ia > 0 ) return(1);
4017  ta += *ta;
4018  tb += *tb;
4019  }
4020  if ( type == THETA ) return(1);
4021  else return(0); /* THETA2 */
4022 }
4023 
4024 /*
4025  #] DoTheta :
4026  #[ DoDelta :
4027 */
4028 
4029 WORD DoDelta(WORD *t)
4030 {
4031  WORD k, *r1, *r2, *tstop, isnum, isnum2, type = *t;
4032  if ( AC.BracketNormalize ) return(-1);
4033  k = t[1];
4034  if ( k <= FUNHEAD ) goto argzero;
4035  if ( k == FUNHEAD+ARGHEAD && t[FUNHEAD] == ARGHEAD ) goto argzero;
4036  tstop = t + k;
4037  t += FUNHEAD;
4038  r1 = t;
4039  NEXTARG(r1)
4040  if ( *t < 0 ) {
4041  k = 1;
4042  if ( *t == -SNUMBER ) { isnum = 1; k = t[1]; }
4043  else isnum = 0;
4044  }
4045  else {
4046  k = t[*t-1];
4047  k = ABS(k);
4048  if ( k == *t-ARGHEAD-1 ) isnum = 1;
4049  else isnum = 0;
4050  k = 1;
4051  }
4052  if ( r1 >= tstop ) { /* Single argument */
4053  if ( !isnum ) return(-1);
4054  if ( k == 0 ) goto argzero;
4055  goto argnonzero;
4056  }
4057  r2 = r1;
4058  NEXTARG(r2)
4059  if ( r2 < tstop ) return(-1);
4060  if ( *r1 < 0 ) {
4061  if ( *r1 == -SNUMBER ) { isnum2 = 1; }
4062  else isnum2 = 0;
4063  }
4064  else {
4065  k = r1[*r1-1];
4066  k = ABS(k);
4067  if ( k == *r1-ARGHEAD-1 ) isnum2 = 1;
4068  else isnum2 = 0;
4069  }
4070  if ( isnum != isnum2 ) return(-1);
4071  tstop = r1;
4072  while ( t < tstop && r1 < r2 ) {
4073  if ( *t != *r1 ) {
4074  if ( !isnum ) return(-1);
4075  goto argnonzero;
4076  }
4077  t++; r1++;
4078  }
4079  if ( t != tstop || r1 != r2 ) {
4080  if ( !isnum ) return(-1);
4081  goto argnonzero;
4082  }
4083 argzero:
4084  if ( type == DELTA2 ) return(1);
4085  else return(0);
4086 argnonzero:
4087  if ( type == DELTA2 ) return(0);
4088  else return(1);
4089 }
4090 
4091 /*
4092  #] DoDelta :
4093  #[ DoRevert :
4094 */
4095 
4096 void DoRevert(WORD *fun, WORD *tmp)
4097 {
4098  WORD *t, *r, *m, *to, *tt, *mm, i, j;
4099  to = fun + fun[1];
4100  r = fun + FUNHEAD;
4101  while ( r < to ) {
4102  if ( *r <= 0 ) {
4103  if ( *r == -REVERSEFUNCTION ) {
4104  m = r; mm = m+1;
4105  while ( mm < to ) *m++ = *mm++;
4106  to--;
4107  (fun[1])--;
4108  fun[2] |= DIRTYSYMFLAG;
4109  }
4110  else if ( *r <= -FUNCTION ) r++;
4111  else {
4112  if ( *r == -INDEX && r[1] < MINSPEC ) *r = -VECTOR;
4113  r += 2;
4114  }
4115  }
4116  else {
4117  if ( ( *r > ARGHEAD )
4118  && ( r[ARGHEAD+1] == REVERSEFUNCTION )
4119  && ( *r == (r[ARGHEAD]+ARGHEAD) )
4120  && ( r[ARGHEAD] == (r[ARGHEAD+2]+4) )
4121  && ( *(r+*r-3) == 1 )
4122  && ( *(r+*r-2) == 1 )
4123  && ( *(r+*r-1) == 3 ) ) {
4124  mm = r;
4125  r += ARGHEAD + 1;
4126  tt = r + r[1];
4127  r += FUNHEAD;
4128  m = tmp;
4129  t = r;
4130  j = 0;
4131  while ( t < tt ) {
4132  NEXTARG(t)
4133  j++;
4134  }
4135  while ( --j >= 0 ) {
4136  i = j;
4137  t = r;
4138  while ( --i >= 0 ) {
4139  NEXTARG(t)
4140  }
4141  if ( *t > 0 ) {
4142  i = *t;
4143  NCOPY(m,t,i);
4144  }
4145  else if ( *t <= -FUNCTION ) *m++ = *t++;
4146  else { *m++ = *t++; *m++ = *t++; }
4147  }
4148  i = WORDDIF(m,tmp);
4149  m = tmp;
4150  t = mm;
4151  r = t + *t;
4152  NCOPY(t,m,i);
4153  m = r;
4154  r = t;
4155  i = WORDDIF(to,m);
4156  NCOPY(t,m,i);
4157  fun[1] = WORDDIF(t,fun);
4158  to = t;
4159  fun[2] |= DIRTYSYMFLAG;
4160  }
4161  else r += *r;
4162  }
4163  }
4164 }
4165 
4166 /*
4167  #] DoRevert :
4168  #] Normalize :
4169  #[ DetCommu :
4170 
4171  Determines the number of terms in an expression that contain
4172  noncommuting objects. This can be used to see whether products of
4173  this expression can be evaluated with binomial coefficients.
4174 
4175  We don't try to be fancy. If a term contains noncommuting objects
4176  we are not looking whether they can commute with complete other
4177  terms.
4178 
4179  If the number gets too large we cut it off.
4180 */
4181 
4182 #define MAXNUMBEROFNONCOMTERMS 2
4183 
4184 WORD DetCommu(WORD *terms)
4185 {
4186  WORD *t, *tnext, *tstop;
4187  WORD num = 0;
4188  if ( *terms == 0 ) return(0);
4189  if ( terms[*terms] == 0 ) return(0);
4190  t = terms;
4191  while ( *t ) {
4192  tnext = t + *t;
4193  tstop = tnext - ABS(tnext[-1]);
4194  t++;
4195  while ( t < tstop ) {
4196  if ( *t >= FUNCTION ) {
4197  if ( functions[*t-FUNCTION].commute ) {
4198  num++;
4199  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4200  break;
4201  }
4202  }
4203  else if ( *t == SUBEXPRESSION ) {
4204  if ( cbuf[t[4]].CanCommu[t[2]] ) {
4205  num++;
4206  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4207  break;
4208  }
4209  }
4210  else if ( *t == EXPRESSION ) {
4211  num++;
4212  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4213  break;
4214  }
4215  else if ( *t == DOLLAREXPRESSION ) {
4216 /*
4217  Technically this is not correct. We have to test first
4218  whether this is MODLOCAL (in TFORM) and if so, use the
4219  local version. Anyway, this should be rare to never
4220  occurring because dollars should be replaced.
4221 */
4222  if ( cbuf[AM.dbufnum].CanCommu[t[2]] ) {
4223  num++;
4224  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4225  break;
4226  }
4227  }
4228  t += t[1];
4229  }
4230  t = tnext;
4231  }
4232  return(num);
4233 }
4234 
4235 /*
4236  #] DetCommu :
4237  #[ DoesCommu :
4238 
4239  Determines the number of noncommuting objects in a term.
4240  If the number gets too large we cut it off.
4241 */
4242 
4243 WORD DoesCommu(WORD *term)
4244 {
4245  WORD *tstop;
4246  WORD num = 0;
4247  if ( *term == 0 ) return(0);
4248  tstop = term + *term;
4249  tstop = tstop - ABS(tstop[-1]);
4250  term++;
4251  while ( term < tstop ) {
4252  if ( ( *term >= FUNCTION ) && ( functions[*term-FUNCTION].commute ) ) {
4253  num++;
4254  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4255  }
4256  term += term[1];
4257  }
4258  return(num);
4259 }
4260 
4261 /*
4262  #] DoesCommu :
4263  #[ PolyNormPoly :
4264 
4265  Normalizes a polynomial
4266 */
4267 
4268 #ifdef EVALUATEGCD
4269 WORD *PolyNormPoly (PHEAD WORD *Poly) {
4270 
4271  GETBIDENTITY;
4272  WORD *buffer = AT.WorkPointer;
4273  WORD *p;
4274  if ( NewSort(BHEAD0) ) { Terminate(-1); }
4275  AR.CompareRoutine = (void *)&CompareSymbols;
4276  while ( *Poly ) {
4277  p = Poly + *Poly;
4278  if ( SymbolNormalize(Poly) < 0 ) return(0);
4279  if ( StoreTerm(BHEAD Poly) ) {
4280  AR.CompareRoutine = (void *)&Compare1;
4281  LowerSortLevel();
4282  Terminate(-1);
4283  }
4284  Poly = p;
4285  }
4286  if ( EndSort(BHEAD buffer,1) < 0 ) {
4287  AR.CompareRoutine = (void *)&Compare1;
4288  Terminate(-1);
4289  }
4290  p = buffer;
4291  while ( *p ) p += *p;
4292  AR.CompareRoutine = (void *)&Compare1;
4293  AT.WorkPointer = p + 1;
4294  return(buffer);
4295 }
4296 #endif
4297 
4298 /*
4299  #] PolyNormPoly :
4300  #[ EvaluateGcd :
4301 
4302  Try to evaluate the GCDFUNCTION gcd_.
4303  This function can have a number of arguments which can be numbers
4304  and/or polynomials. If there are objects that aren't SYMBOLS or numbers
4305  it cannot work currently.
4306 
4307  To make this work properly we have to intervene in proces.c
4308  proces.c: if ( Normalize(BHEAD m) ) {
4309 1060 proces.c: if ( Normalize(BHEAD r) ) {
4310 1126?proces.c: if ( Normalize(BHEAD term) ) {
4311  proces.c: if ( Normalize(BHEAD AT.WorkPointer) ) goto PasErr;
4312 2308!proces.c: if ( ( retnorm = Normalize(BHEAD term) ) != 0 ) {
4313  proces.c: ReNumber(BHEAD term); Normalize(BHEAD term);
4314  proces.c: if ( Normalize(BHEAD v) ) Terminate(-1);
4315  proces.c: if ( Normalize(BHEAD w) ) { LowerSortLevel(); goto PolyCall; }
4316  proces.c: if ( Normalize(BHEAD term) ) goto PolyCall;
4317 */
4318 #ifdef EVALUATEGCD
4319 
4320 WORD *EvaluateGcd(PHEAD WORD *subterm)
4321 {
4322  GETBIDENTITY
4323  WORD *oldworkpointer = AT.WorkPointer, *work1, *work2, *work3;
4324  WORD *t, *tt, *ttt, *t1, *t2, *t3, *t4, *tstop;
4325  WORD ct, nnum;
4326  UWORD gcdnum, stor;
4327  WORD *lnum=AT.n_llnum+1;
4328  WORD *num1, *num2, *num3, *den1, *den2, *den3;
4329  WORD sizenum1, sizenum2, sizenum3, sizeden1, sizeden2, sizeden3;
4330  int i, isnumeric = 0, numarg = 0 /*, sizearg */;
4331  LONG size;
4332 /*
4333  Step 1: Look for -SNUMBER or -SYMBOL arguments.
4334  If encountered, treat everybody with it.
4335 */
4336  tt = subterm + subterm[1]; t = subterm + FUNHEAD;
4337 
4338  while ( t < tt ) {
4339  numarg++;
4340  if ( *t == -SNUMBER ) {
4341  if ( t[1] == 0 ) {
4342 gcdzero:;
4343  MLOCK(ErrorMessageLock);
4344  MesPrint("Trying to take the GCD involving a zero term.");
4345  MUNLOCK(ErrorMessageLock);
4346  return(0);
4347  }
4348  gcdnum = ABS(t[1]);
4349  t1 = subterm + FUNHEAD;
4350  while ( gcdnum > 1 && t1 < tt ) {
4351  if ( *t1 == -SNUMBER ) {
4352  stor = ABS(t1[1]);
4353  if ( stor == 0 ) goto gcdzero;
4354  if ( GcdLong(BHEAD (UWORD *)&stor,1,(UWORD *)&gcdnum,1,
4355  (UWORD *)lnum,&nnum) ) goto FromGCD;
4356  gcdnum = lnum[0];
4357  t1 += 2;
4358  continue;
4359  }
4360  else if ( *t1 == -SYMBOL ) goto gcdisone;
4361  else if ( *t1 < 0 ) goto gcdillegal;
4362 /*
4363  Now we have to go through all the terms in the argument.
4364  This includes long numbers.
4365 */
4366  ttt = t1 + *t1;
4367  ct = *ttt; *ttt = 0;
4368  if ( t1[1] != 0 ) { /* First normalize the argument */
4369  t1 = PolyNormPoly(BHEAD t1+ARGHEAD);
4370  }
4371  else t1 += ARGHEAD;
4372  while ( *t1 ) {
4373  t1 += *t1;
4374  i = ABS(t1[-1]);
4375  t2 = t1 - i;
4376  i = (i-1)/2;
4377  t3 = t2+i-1;
4378  while ( t3 > t2 && *t3 == 0 ) { t3--; i--; }
4379  if ( GcdLong(BHEAD (UWORD *)t2,(WORD)i,(UWORD *)&gcdnum,1,
4380  (UWORD *)lnum,&nnum) ) {
4381  *ttt = ct;
4382  goto FromGCD;
4383  }
4384  gcdnum = lnum[0];
4385  if ( gcdnum == 1 ) {
4386  *ttt = ct;
4387  goto gcdisone;
4388  }
4389  }
4390  *ttt = ct;
4391  t1 = ttt;
4392  AT.WorkPointer = oldworkpointer;
4393  }
4394  if ( gcdnum == 1 ) goto gcdisone;
4395  oldworkpointer[0] = 4;
4396  oldworkpointer[1] = gcdnum;
4397  oldworkpointer[2] = 1;
4398  oldworkpointer[3] = 3;
4399  oldworkpointer[4] = 0;
4400  AT.WorkPointer = oldworkpointer + 5;
4401  return(oldworkpointer);
4402  }
4403  else if ( *t == -SYMBOL ) {
4404  t1 = subterm + FUNHEAD;
4405  i = t[1];
4406  while ( t1 < tt ) {
4407  if ( *t1 == -SNUMBER ) goto gcdisone;
4408  if ( *t1 == -SYMBOL ) {
4409  if ( t1[1] != i ) goto gcdisone;
4410  t1 += 2; continue;
4411  }
4412  if ( *t1 < 0 ) goto gcdillegal;
4413  ttt = t1 + *t1;
4414  ct = *ttt; *ttt = 0;
4415  if ( t1[1] != 0 ) { /* First normalize the argument */
4416  t2 = PolyNormPoly(BHEAD t1+ARGHEAD);
4417  }
4418  else t2 = t1 + ARGHEAD;
4419  while ( *t2 ) {
4420  t3 = t2+1;
4421  t2 = t2 + *t2;
4422  tstop = t2 - ABS(t2[-1]);
4423  while ( t3 < tstop ) {
4424  if ( *t3 != SYMBOL ) {
4425  *ttt = ct;
4426  goto gcdillegal;
4427  }
4428  t4 = t3 + 2;
4429  t3 += t3[1];
4430  while ( t4 < t3 ) {
4431  if ( *t4 == i && t4[1] > 0 ) goto nextterminarg;
4432  t4 += 2;
4433  }
4434  }
4435  *ttt = ct;
4436  goto gcdisone;
4437 nextterminarg:;
4438  }
4439  *ttt = ct;
4440  t1 = ttt;
4441  AT.WorkPointer = oldworkpointer;
4442  }
4443  oldworkpointer[0] = 8;
4444  oldworkpointer[1] = SYMBOL;
4445  oldworkpointer[2] = 4;
4446  oldworkpointer[3] = t[1];
4447  oldworkpointer[4] = 1;
4448  oldworkpointer[5] = 1;
4449  oldworkpointer[6] = 1;
4450  oldworkpointer[7] = 3;
4451  oldworkpointer[8] = 0;
4452  AT.WorkPointer = oldworkpointer+9;
4453  return(oldworkpointer);
4454  }
4455  else if ( *t < 0 ) {
4456 gcdillegal:;
4457  MLOCK(ErrorMessageLock);
4458  MesPrint("Illegal object in gcd_ function. Object not a number or a symbol.");
4459  MUNLOCK(ErrorMessageLock);
4460  goto FromGCD;
4461  }
4462  else if ( ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4463  else if ( t[1] != 0 ) {
4464  ttt = t + *t; ct = *ttt; *ttt = 0;
4465  t = PolyNormPoly(BHEAD t+ARGHEAD);
4466  *ttt = ct;
4467  if ( t[*t] == 0 && ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4468  AT.WorkPointer = oldworkpointer;
4469  t = ttt;
4470  }
4471  t += *t;
4472  }
4473 /*
4474  At this point there are only generic arguments.
4475  There are however still two cases:
4476  1: There is an argument that is purely numerical
4477  In that case we have to take the gcd of all coefficients
4478  2: All arguments are nontrivial polynomials.
4479  Here we don't worry so much about the factor. (???)
4480  We know whether case 1 occurs when isnumeric > 0.
4481  We can look up numarg to get a good starting value.
4482 */
4483  AT.WorkPointer = oldworkpointer;
4484  if ( isnumeric ) {
4485  t = subterm + FUNHEAD;
4486  for ( i = 1; i < isnumeric; i++ ) {
4487  NEXTARG(t);
4488  }
4489  if ( t[1] != 0 ) { /* First normalize the argument */
4490  ttt = t + *t; ct = *ttt; *ttt = 0;
4491  t = PolyNormPoly(BHEAD t+ARGHEAD);
4492  *ttt = ct;
4493  }
4494  t += *t;
4495  i = (ABS(t[-1])-1)/2;
4496  den1 = t - 1 - i;
4497  num1 = den1 - i;
4498  sizenum1 = sizeden1 = i;
4499  while ( sizenum1 > 1 && num1[sizenum1-1] == 0 ) sizenum1--;
4500  while ( sizeden1 > 1 && den1[sizeden1-1] == 0 ) sizeden1--;
4501  work1 = AT.WorkPointer+1; work2 = work1+sizenum1;
4502  for ( i = 0; i < sizenum1; i++ ) work1[i] = num1[i];
4503  for ( i = 0; i < sizeden1; i++ ) work2[i] = den1[i];
4504  num1 = work1; den1 = work2;
4505  AT.WorkPointer = work2 = work2 + sizeden1;
4506  t = subterm + FUNHEAD;
4507  while ( t < tt ) {
4508  ttt = t + *t; ct = *ttt; *ttt = 0;
4509  if ( t[1] != 0 ) {
4510  t = PolyNormPoly(BHEAD t+ARGHEAD);
4511  }
4512  else t += ARGHEAD;
4513  while ( *t ) {
4514  t += *t;
4515  i = (ABS(t[-1])-1)/2;
4516  den2 = t - 1 - i;
4517  num2 = den2 - i;
4518  sizenum2 = sizeden2 = i;
4519  while ( sizenum2 > 1 && num2[sizenum2-1] == 0 ) sizenum2--;
4520  while ( sizeden2 > 1 && den2[sizeden2-1] == 0 ) sizeden2--;
4521  num3 = AT.WorkPointer;
4522  if ( GcdLong(BHEAD (UWORD *)num2,sizenum2,(UWORD *)num1,sizenum1,
4523  (UWORD *)num3,&sizenum3) ) goto FromGCD;
4524  sizenum1 = sizenum3;
4525  for ( i = 0; i < sizenum1; i++ ) num1[i] = num3[i];
4526  den3 = AT.WorkPointer;
4527  if ( GcdLong(BHEAD (UWORD *)den2,sizeden2,(UWORD *)den1,sizeden1,
4528  (UWORD *)den3,&sizeden3) ) goto FromGCD;
4529  sizeden1 = sizeden3;
4530  for ( i = 0; i < sizeden1; i++ ) den1[i] = den3[i];
4531  if ( sizenum1 == 1 && num1[0] == 1 && sizeden1 == 1 && den1[1] == 1 )
4532  goto gcdisone;
4533  }
4534  *ttt = ct;
4535  t = ttt;
4536  AT.WorkPointer = work2;
4537  }
4538  AT.WorkPointer = oldworkpointer;
4539 /*
4540  Now copy the GCD to the 'output'
4541 */
4542  if ( sizenum1 > sizeden1 ) {
4543  while ( sizenum1 > sizeden1 ) den1[sizeden1++] = 0;
4544  }
4545  else if ( sizenum1 < sizeden1 ) {
4546  while ( sizenum1 < sizeden1 ) num1[sizenum1++] = 0;
4547  }
4548  t = oldworkpointer;
4549  i = 2*sizenum1+1;
4550  *t++ = i+1;
4551  if ( num1 != t ) { NCOPY(t,num1,sizenum1); }
4552  else t += sizenum1;
4553  if ( den1 != t ) { NCOPY(t,den1,sizeden1); }
4554  else t += sizeden1;
4555  *t++ = i;
4556  *t++ = 0;
4557  AT.WorkPointer = t;
4558  return(oldworkpointer);
4559  }
4560 /*
4561  Now the real stuff with only polynomials.
4562  Pick up the shortest term to start.
4563  We are a bit brutish about this.
4564 */
4565  t = subterm + FUNHEAD;
4566  AT.WorkPointer += AM.MaxTer/sizeof(WORD);
4567  work2 = AT.WorkPointer;
4568 /*
4569  sizearg = subterm[1];
4570 */
4571  i = 0; work3 = 0;
4572  while ( t < tt ) {
4573  i++;
4574  work1 = AT.WorkPointer;
4575  ttt = t + *t; ct = *ttt; *ttt = 0;
4576  t = PolyNormPoly(BHEAD t+ARGHEAD);
4577  if ( *work1 < AT.WorkPointer-work1 ) {
4578 /*
4579  sizearg = AT.WorkPointer-work1;
4580 */
4581  numarg = i;
4582  work3 = work1;
4583  }
4584  *ttt = ct; t = ttt;
4585  }
4586  *AT.WorkPointer++ = 0;
4587 /*
4588  We have properly normalized arguments and the shortest is indicated in work3
4589 */
4590  work1 = work3;
4591  while ( *work2 ) {
4592  if ( work2 != work3 ) {
4593  work1 = PolyGCD2(BHEAD work1,work2);
4594  }
4595  while ( *work2 ) work2 += *work2;
4596  work2++;
4597  }
4598  work2 = work1;
4599  while ( *work2 ) work2 += *work2;
4600  size = work2 - work1 + 1;
4601  t = oldworkpointer;
4602  NCOPY(t,work1,size);
4603  AT.WorkPointer = t;
4604  return(oldworkpointer);
4605 
4606 gcdisone:;
4607  oldworkpointer[0] = 4;
4608  oldworkpointer[1] = 1;
4609  oldworkpointer[2] = 1;
4610  oldworkpointer[3] = 3;
4611  oldworkpointer[4] = 0;
4612  AT.WorkPointer = oldworkpointer+5;
4613  return(oldworkpointer);
4614 FromGCD:
4615  MLOCK(ErrorMessageLock);
4616  MesCall("EvaluateGcd");
4617  MUNLOCK(ErrorMessageLock);
4618  return(0);
4619 }
4620 
4621 #endif
4622 
4623 /*
4624  #] EvaluateGcd :
4625  #[ DropCoefficient :
4626 */
4627 
4628 void DropCoefficient(PHEAD WORD *term)
4629 {
4630  GETBIDENTITY
4631  WORD *t = term + *term;
4632  WORD n, na;
4633  n = t[-1]; na = ABS(n);
4634  t -= na;
4635  if ( n == 3 && t[0] == 1 && t[1] == 1 ) return;
4636  *AN.RepPoint = 1;
4637  t[0] = 1; t[1] = 1; t[2] = 3;
4638  *term -= (na-3);
4639 }
4640 
4641 /*
4642  #] DropCoefficient :
4643  #[ DropSymbols :
4644 */
4645 
4646 void DropSymbols(PHEAD WORD *term)
4647 {
4648  GETBIDENTITY
4649  WORD *tend = term + *term, *t1, *t2, *tstop;
4650  tstop = tend - ABS(tend[-1]);
4651  t1 = term+1;
4652  while ( t1 < tstop ) {
4653  if ( *t1 == SYMBOL ) {
4654  *AN.RepPoint = 1;
4655  t2 = t1+t1[1];
4656  while ( t2 < tend ) *t1++ = *t2++;
4657  *term = t1 - term;
4658  break;
4659  }
4660  t1 += t1[1];
4661  }
4662 }
4663 
4664 /*
4665  #] DropSymbols :
4666  #[ SymbolNormalize :
4667 */
4676 int SymbolNormalize(WORD *term)
4677 {
4678  WORD buffer[7*NORMSIZE], *t, *b, *bb, *tt, *m, *tstop;
4679  int i;
4680  b = buffer;
4681  *b++ = SYMBOL; *b++ = 2;
4682  t = term + *term;
4683  tstop = t - ABS(t[-1]);
4684  t = term + 1;
4685  while ( t < tstop ) { /* Step 1: collect symbols */
4686  if ( *t == SYMBOL && t < tstop ) {
4687  for ( i = 2; i < t[1]; i += 2 ) {
4688  bb = buffer+2;
4689  while ( bb < b ) {
4690  if ( bb[0] == t[i] ) { /* add powers */
4691  bb[1] += t[i+1];
4692  if ( bb[1] > MAXPOWER || bb[1] < -MAXPOWER ) {
4693  MLOCK(ErrorMessageLock);
4694  MesPrint("Power in SymbolNormalize out of range");
4695  MUNLOCK(ErrorMessageLock);
4696  return(-1);
4697  }
4698  if ( bb[1] == 0 ) {
4699  b -= 2;
4700  while ( bb < b ) {
4701  bb[0] = bb[2]; bb[1] = bb[3]; bb += 2;
4702  }
4703  }
4704  goto Nexti;
4705  }
4706  else if ( bb[0] > t[i] ) { /* insert it */
4707  m = b;
4708  while ( m > bb ) { m[1] = m[-1]; m[0] = m[-2]; m -= 2; }
4709  b += 2;
4710  bb[0] = t[i];
4711  bb[1] = t[i+1];
4712  goto Nexti;
4713  }
4714  bb += 2;
4715  }
4716  if ( bb >= b ) { /* add it to the end */
4717  *b++ = t[i]; *b++ = t[i+1];
4718  }
4719 Nexti:;
4720  }
4721  }
4722  else {
4723  MLOCK(ErrorMessageLock);
4724  MesPrint("Illegal term in SymbolNormalize");
4725  MUNLOCK(ErrorMessageLock);
4726  return(-1);
4727  }
4728  t += t[1];
4729  }
4730  buffer[1] = b - buffer;
4731 /*
4732  Veto negative powers
4733 */
4734  b = buffer; bb = b + b[1]; b += 3;
4735  while ( b < bb ) {
4736  if ( *b < 0 ) {
4737  MLOCK(ErrorMessageLock);
4738  MesPrint("Negative power in SymbolNormalize");
4739  MUNLOCK(ErrorMessageLock);
4740  return(-1);
4741  }
4742  b += 2;
4743  }
4744 /*
4745  Now we use the fact that the new term will not be longer than the old one
4746  Actually it should be shorter when there is more than one subterm!
4747  Copy back.
4748 */
4749  i = buffer[1];
4750  b = buffer; tt = term + 1;
4751  if ( i > 2 ) { NCOPY(tt,b,i) }
4752  if ( tt < tstop ) {
4753  i = term[*term-1];
4754  if ( i < 0 ) i = -i;
4755  *term -= (tstop-tt);
4756  NCOPY(tt,tstop,i)
4757  }
4758  return(0);
4759 }
4760 
4761 /*
4762  #] SymbolNormalize :
4763 */
#define PHEAD
Definition: ftypes.h:56
int SymbolNormalize(WORD *term)
Definition: normal.c:4676
Definition: structs.h:908
WORD * Pointer
Definition: structs.h:911
WORD StoreTerm(PHEAD WORD *)
Definition: sort.c:4070
WORD ** rhs
Definition: structs.h:913
WORD Compare1(PHEAD WORD *, WORD *, WORD)
Definition: sort.c:2397
VOID LowerSortLevel()
Definition: sort.c:4435
WORD * Buffer
Definition: structs.h:909
WORD NewSort(PHEAD0)
Definition: sort.c:553
WORD NextPrime(PHEAD WORD)
Definition: reken.c:3629
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:2865
WORD * Top
Definition: structs.h:910
int CompareSymbols(PHEAD WORD *, WORD *, WORD)
Definition: sort.c:2818
WORD CompCoef(WORD *, WORD *)
Definition: reken.c:3012
LONG EndSort(PHEAD WORD *, int)
Definition: sort.c:632