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