FORM  4.3
normal.c
Go to the documentation of this file.
1 
10 /* #[ License : */
11 /*
12  * Copyright (C) 1984-2022 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 MOEBIUS:
1963 /*
1964  Numerical function giving -1,0,1 or no evaluation
1965 */
1966  if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER
1967  && t[FUNHEAD+1] > 0 ) {
1968  WORD val = Moebius(BHEAD t[FUNHEAD+1]);
1969  if ( val == 0 ) goto NormZero;
1970  if ( val < 0 ) ncoef = -ncoef;
1971  }
1972  else {
1973  pcom[ncom++] = t;
1974  }
1975  break;
1976  case LNUMBER :
1977  ncoef = REDLENG(ncoef);
1978  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(t+3),t[2]) ) goto FromNorm;
1979  ncoef = INCLENG(ncoef);
1980  break;
1981  case SNUMBER :
1982  if ( t[2] < 0 ) {
1983  t[2] = -t[2];
1984  if ( t[3] & 1 ) ncoef = -ncoef;
1985  }
1986  else if ( t[2] == 0 ) {
1987  if ( t[3] < 0 ) goto NormInf;
1988  goto NormZero;
1989  }
1990  lnum[0] = t[2];
1991  nnum = 1;
1992  if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) ) goto FromNorm;
1993  ncoef = REDLENG(ncoef);
1994  if ( t[3] < 0 ) {
1995  if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1996  goto FromNorm;
1997  }
1998  else if ( t[3] > 0 ) {
1999  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2000  goto FromNorm;
2001  }
2002  ncoef = INCLENG(ncoef);
2003  break;
2004  case GAMMA :
2005  case GAMMAI :
2006  case GAMMAFIVE :
2007  case GAMMASIX :
2008  case GAMMASEVEN :
2009  if ( t[1] == FUNHEAD ) {
2010  MLOCK(ErrorMessageLock);
2011  MesPrint("Gamma matrix without spin line encountered.");
2012  MUNLOCK(ErrorMessageLock);
2013  goto NormMin;
2014  }
2015  pnco[nnco++] = t;
2016  t += FUNHEAD+1;
2017  goto ScanCont;
2018  case LEVICIVITA :
2019  peps[neps++] = t;
2020  if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG ) {
2021  t[2] &= ~DIRTYFLAG;
2022  t[2] |= DIRTYSYMFLAG;
2023  }
2024  t += FUNHEAD;
2025 ScanCont: while ( t < r ) {
2026  if ( *t >= AM.OffsetIndex &&
2027  ( *t >= AM.DumInd || ( *t < AM.WilInd &&
2028  indices[*t-AM.OffsetIndex].dimension ) ) )
2029  pcon[ncon++] = t;
2030  t++;
2031  }
2032  break;
2033  case EXPONENT :
2034  {
2035  WORD *rr;
2036  k = 1;
2037  rr = t + FUNHEAD;
2038  if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) )
2039  k = 0;
2040  if ( *rr == -SNUMBER && rr[1] == 1 ) break;
2041  if ( *rr <= -FUNCTION ) k = *rr;
2042  NEXTARG(rr)
2043  if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) ) {
2044  if ( k == 0 ) goto NormZZ;
2045  break;
2046  }
2047  if ( *rr == -SNUMBER && rr[1] > 0 && rr[1] < MAXPOWER && k < 0 ) {
2048  k = -k;
2049  if ( functions[k-FUNCTION].commute ) {
2050  for ( i = 0; i < rr[1]; i++ ) pnco[nnco++] = rr-1;
2051  }
2052  else {
2053  for ( i = 0; i < rr[1]; i++ ) pcom[ncom++] = rr-1;
2054  }
2055  break;
2056  }
2057  if ( k == 0 ) goto NormZero;
2058  if ( t[FUNHEAD] == -SYMBOL && *rr == -SNUMBER && t[1] == FUNHEAD+4 ) {
2059  if ( rr[1] < MAXPOWER ) {
2060  t[FUNHEAD+2] = t[FUNHEAD+1]; t += FUNHEAD+2;
2061  from = m;
2062  goto NextSymbol;
2063  }
2064  }
2065 /*
2066  if ( ( t[FUNHEAD] > 0 && t[FUNHEAD+1] != 0 )
2067  || ( *rr > 0 && rr[1] != 0 ) ) {}
2068  else
2069 */
2070  t[2] &= ~DIRTYSYMFLAG;
2071 
2072  pnco[nnco++] = t;
2073  }
2074  break;
2075  case DENOMINATOR :
2076  t[2] &= ~DIRTYSYMFLAG;
2077  pden[nden++] = t;
2078  pnco[nnco++] = t;
2079  break;
2080  case INDEX :
2081  t += 2;
2082  do {
2083  if ( *t == 0 || *t == AM.vectorzero ) goto NormZero;
2084  if ( *t > 0 && *t < AM.OffsetIndex ) {
2085  lnum[0] = *t++;
2086  nnum = 1;
2087  ncoef = REDLENG(ncoef);
2088  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2089  goto FromNorm;
2090  ncoef = INCLENG(ncoef);
2091  }
2092  else if ( *t == NOINDEX ) t++;
2093  else pind[nind++] = *t++;
2094  } while ( t < r );
2095  break;
2096  case SUBEXPRESSION :
2097  if ( t[3] == 0 ) break;
2098  case EXPRESSION :
2099  goto RegEnd;
2100  case ROOTFUNCTION :
2101 /*
2102  Tries to take the n-th root inside the rationals
2103  If this is not possible, it clears all flags and
2104  hence tries no more.
2105  Notation:
2106  root_(power(=integer),(rational)number)
2107 */
2108  { WORD nc;
2109  if ( t[2] == 0 ) goto defaultcase;
2110  if ( t[FUNHEAD] != -SNUMBER || t[FUNHEAD+1] < 0 ) goto defaultcase;
2111  if ( t[FUNHEAD+2] == -SNUMBER ) {
2112  if ( t[FUNHEAD+1] == 0 && t[FUNHEAD+3] == 0 ) goto NormZZ;
2113  if ( t[FUNHEAD+1] == 0 ) break;
2114  if ( t[FUNHEAD+3] < 0 ) {
2115  AT.WorkPointer[0] = -t[FUNHEAD+3];
2116  nc = -1;
2117  }
2118  else {
2119  AT.WorkPointer[0] = t[FUNHEAD+3];
2120  nc = 1;
2121  }
2122  AT.WorkPointer[1] = 1;
2123  }
2124  else if ( t[FUNHEAD+2] == t[1]-FUNHEAD-2
2125  && t[FUNHEAD+2] == t[FUNHEAD+2+ARGHEAD]+ARGHEAD
2126  && ABS(t[t[1]-1]) == t[FUNHEAD+2+ARGHEAD] - 1 ) {
2127  WORD *r1, *r2;
2128  if ( t[FUNHEAD+1] == 0 ) break;
2129  i = t[t[1]-1]; r1 = t + FUNHEAD+ARGHEAD+3;
2130  nc = REDLENG(i);
2131  i = ABS(i) - 1;
2132  r2 = AT.WorkPointer;
2133  while ( --i >= 0 ) *r2++ = *r1++;
2134  }
2135  else goto defaultcase;
2136  if ( TakeRatRoot((UWORD *)AT.WorkPointer,&nc,t[FUNHEAD+1]) ) {
2137  t[2] = 0;
2138  goto defaultcase;
2139  }
2140  ncoef = REDLENG(ncoef);
2141  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,nc) )
2142  goto FromNorm;
2143  if ( nc < 0 ) nc = -nc;
2144  if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(AT.WorkPointer+nc),nc) )
2145  goto FromNorm;
2146  ncoef = INCLENG(ncoef);
2147  }
2148  break;
2149  case RANDOMFUNCTION :
2150  {
2151  WORD nnc, nc, nca, nr;
2152  UWORD xx;
2153 /*
2154  Needs one positive integer argument.
2155  returns (wranf()%argument)+1.
2156  We may call wranf several times to paste UWORDS together
2157  when we need long numbers.
2158  We make little errors when taking the % operator
2159  (not 100% uniform). We correct for that by redoing the
2160  calculation in the (unlikely) case that we are in leftover area
2161 */
2162  if ( t[1] == FUNHEAD ) goto defaultcase;
2163  if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER &&
2164  t[FUNHEAD+1] > 0 ) {
2165  if ( t[FUNHEAD+1] == 1 ) break;
2166 redoshort:
2167  ((UWORD *)AT.WorkPointer)[0] = wranf(BHEAD0);
2168  ((UWORD *)AT.WorkPointer)[1] = wranf(BHEAD0);
2169  nr = 2;
2170  if ( ((UWORD *)AT.WorkPointer)[1] == 0 ) {
2171  nr = 1;
2172  if ( ((UWORD *)AT.WorkPointer)[0] == 0 ) {
2173  nr = 0;
2174  }
2175  }
2176  xx = (UWORD)(t[FUNHEAD+1]);
2177  if ( nr ) {
2178  DivLong((UWORD *)AT.WorkPointer,nr
2179  ,&xx,1
2180  ,((UWORD *)AT.WorkPointer)+4,&nnc
2181  ,((UWORD *)AT.WorkPointer)+2,&nc);
2182  ((UWORD *)AT.WorkPointer)[4] = 0;
2183  ((UWORD *)AT.WorkPointer)[5] = 0;
2184  ((UWORD *)AT.WorkPointer)[6] = 1;
2185  DivLong((UWORD *)AT.WorkPointer+4,3
2186  ,&xx,1
2187  ,((UWORD *)AT.WorkPointer)+9,&nnc
2188  ,((UWORD *)AT.WorkPointer)+7,&nca);
2189  AddLong((UWORD *)AT.WorkPointer+4,3
2190  ,((UWORD *)AT.WorkPointer)+7,-nca
2191  ,((UWORD *)AT.WorkPointer)+9,&nnc);
2192  if ( BigLong((UWORD *)AT.WorkPointer,nr
2193  ,((UWORD *)AT.WorkPointer)+9,nnc) >= 0 ) goto redoshort;
2194  }
2195  else nc = 0;
2196  if ( nc == 0 ) {
2197  AT.WorkPointer[2] = (WORD)xx;
2198  nc = 1;
2199  }
2200  ncoef = REDLENG(ncoef);
2201  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+2,nc) )
2202  goto FromNorm;
2203  ncoef = INCLENG(ncoef);
2204  }
2205  else if ( t[FUNHEAD] > 0 && t[1] == t[FUNHEAD]+FUNHEAD
2206  && ABS(t[t[1]-1]) == t[FUNHEAD]-1-ARGHEAD && t[t[1]-1] > 0 ) {
2207  WORD nna, nnb, nni, nnb2, nnb2a;
2208  UWORD *nnt;
2209  nna = t[t[1]-1];
2210  nnb2 = nna-1;
2211  nnb = nnb2/2;
2212  nnt = (UWORD *)(t+t[1]-1-nnb); /* start of denominator */
2213  if ( *nnt != 1 ) goto defaultcase;
2214  for ( nni = 1; nni < nnb; nni++ ) {
2215  if ( nnt[nni] != 0 ) goto defaultcase;
2216  }
2217  nnt = (UWORD *)(t + FUNHEAD + ARGHEAD + 1);
2218 
2219  for ( nni = 0; nni < nnb2; nni++ ) {
2220  ((UWORD *)AT.WorkPointer)[nni] = wranf(BHEAD0);
2221  }
2222  nnb2a = nnb2;
2223  while ( nnb2a > 0 && ((UWORD *)AT.WorkPointer)[nnb2a-1] == 0 ) nnb2a--;
2224  if ( nnb2a > 0 ) {
2225  DivLong((UWORD *)AT.WorkPointer,nnb2a
2226  ,nnt,nnb
2227  ,((UWORD *)AT.WorkPointer)+2*nnb2,&nnc
2228  ,((UWORD *)AT.WorkPointer)+nnb2,&nc);
2229  for ( nni = 0; nni < nnb2; nni++ ) {
2230  ((UWORD *)AT.WorkPointer)[nni+2*nnb2] = 0;
2231  }
2232  ((UWORD *)AT.WorkPointer)[3*nnb2] = 1;
2233  DivLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2234  ,nnt,nnb
2235  ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc
2236  ,((UWORD *)AT.WorkPointer)+3*nnb2+1,&nca);
2237  AddLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2238  ,((UWORD *)AT.WorkPointer)+3*nnb2+1,-nca
2239  ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc);
2240  if ( BigLong((UWORD *)AT.WorkPointer,nnb2a
2241  ,((UWORD *)AT.WorkPointer)+4*nnb2+1,nnc) >= 0 ) goto redoshort;
2242  }
2243  else nc = 0;
2244  if ( nc == 0 ) {
2245  for ( nni = 0; nni < nnb; nni++ ) {
2246  ((UWORD *)AT.WorkPointer)[nnb2+nni] = nnt[nni];
2247  }
2248  nc = nnb;
2249  }
2250  ncoef = REDLENG(ncoef);
2251  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+nnb2,nc) )
2252  goto FromNorm;
2253  ncoef = INCLENG(ncoef);
2254  }
2255  else goto defaultcase;
2256  }
2257  break;
2258  case RANPERM :
2259  if ( *t == RANPERM && t[1] > FUNHEAD && t[FUNHEAD] <= -FUNCTION ) {
2260  WORD **pwork;
2261  WORD *mm, *ww, *ow = AT.WorkPointer;
2262  WORD *Array, *targ, *argstop, narg = 0, itot;
2263  int ie;
2264  argstop = t+t[1];
2265  targ = t+FUNHEAD+1;
2266  while ( targ < argstop ) {
2267  narg++; NEXTARG(targ);
2268  }
2269  WantAddPointers(narg);
2270  pwork = AT.pWorkSpace + AT.pWorkPointer;
2271  targ = t+FUNHEAD+1; narg = 0;
2272  while ( targ < argstop ) {
2273  pwork[narg++] = targ;
2274  NEXTARG(targ);
2275  }
2276 /*
2277  Make a random permutation of the numbers 0,...,narg-1
2278  The following code works also for narg == 0 and narg == 1
2279 */
2280  ow = AT.WorkPointer;
2281  Array = AT.WorkPointer;
2282  AT.WorkPointer += narg;
2283  for ( i = 0; i < narg; i++ ) Array[i] = i;
2284  for ( i = 2; i <= narg; i++ ) {
2285  itot = (WORD)(iranf(BHEAD i));
2286  for ( j = 0; j < itot; j++ ) CYCLE1(WORD,Array,i)
2287  }
2288  mm = AT.WorkPointer;
2289  *mm++ = -t[FUNHEAD];
2290  *mm++ = t[1] - 1;
2291  for ( ie = 2; ie < FUNHEAD; ie++ ) *mm++ = t[ie];
2292  for ( i = 0; i < narg; i++ ) {
2293  ww = pwork[Array[i]];
2294  CopyArg(mm,ww);
2295  }
2296  mm = AT.WorkPointer; t++; ww = t;
2297  i = mm[1]; NCOPY(ww,mm,i)
2298  AT.WorkPointer = ow;
2299  goto TryAgain;
2300  }
2301  pnco[nnco++] = t;
2302  break;
2303  case PUTFIRST : /* First argument should be a function, second a number */
2304  if ( ( t[2] & DIRTYFLAG ) != 0 && t[FUNHEAD] <= -FUNCTION
2305  && t[FUNHEAD+1] == -SNUMBER && t[FUNHEAD+2] > 0 ) {
2306  WORD *rr = t+t[1], *mm = t+FUNHEAD+3, *tt, *tt1, *tt2, num = 0;
2307 /*
2308  now count the arguments. If not enough: no action.
2309 */
2310  while ( mm < rr ) { num++; NEXTARG(mm); }
2311  if ( num < t[FUNHEAD+2] ) { pnco[nnco++] = t; break; }
2312  *t = -t[FUNHEAD]; mm = t+FUNHEAD+3;
2313  i = t[FUNHEAD+2];
2314  while ( --i > 0 ) { NEXTARG(mm); }
2315  tt = TermMalloc("Select_"); /* Move selected out of the way */
2316  tt1 = tt;
2317  if ( *mm > 0 ) {
2318  for ( i = 0; i < *mm; i++ ) *tt1++ = mm[i];
2319  }
2320  else if ( *mm <= -FUNCTION ) { *tt1++ = *mm; }
2321  else { *tt1++ = mm[0]; *tt1++ = mm[1]; }
2322  tt2 = t+FUNHEAD+3;
2323  while ( tt2 < mm ) *tt1++ = *tt2++;
2324  i = tt1-tt; tt1 = tt; tt2 = t+FUNHEAD;
2325  NCOPY(tt2,tt1,i);
2326  TermFree(tt,"Select_");
2327  NEXTARG(mm);
2328  while ( mm < rr ) *tt2++ = *mm++;
2329  t[1] = tt2 - t;
2330  rr = term + *term;
2331  while ( mm < rr ) *tt2++ = *mm++;
2332  *term = tt2-term;
2333  goto Restart;
2334  }
2335  else pnco[nnco++] = t;
2336  break;
2337  case INTFUNCTION :
2338 /*
2339  Can be resolved if the first argument is a number
2340  and the second argument either doesn't exist or has
2341  the value +1, 0, -1
2342  +1 : rounding up
2343  0 : rounding towards zero
2344  -1 : rounding down (same as no argument)
2345 */
2346  if ( t[1] <= FUNHEAD ) break;
2347  {
2348  WORD *rr, den, num;
2349  to = t + FUNHEAD;
2350  if ( *to > 0 ) {
2351  if ( *to == ARGHEAD ) goto NormZero;
2352  rr = to + *to;
2353  i = rr[-1];
2354  j = ABS(i);
2355  if ( to[ARGHEAD] != j+1 ) goto NoInteg;
2356  if ( rr >= r ) k = -1;
2357  else if ( *rr == ARGHEAD ) { k = 0; rr += ARGHEAD; }
2358  else if ( *rr == -SNUMBER ) { k = rr[1]; rr += 2; }
2359  else goto NoInteg;
2360  if ( rr != r ) goto NoInteg;
2361  if ( k > 1 || k < -1 ) goto NoInteg;
2362  to += ARGHEAD+1;
2363  j = (j-1) >> 1;
2364  i = ( i < 0 ) ? -j: j;
2365  UnPack((UWORD *)to,i,&den,&num);
2366 /*
2367  Potentially the use of NoScrat2 is unsafe.
2368  It makes the routine not reentrant, but because it is
2369  used only locally and because we only call the
2370  low level routines DivLong and AddLong which never
2371  make calls involving Normalize, things are OK after all
2372 */
2373  if ( AN.NoScrat2 == 0 ) {
2374  AN.NoScrat2 = (UWORD *)Malloc1((AM.MaxTal+2)*sizeof(UWORD),"Normalize");
2375  }
2376  if ( DivLong((UWORD *)to,num,(UWORD *)(to+j),den
2377  ,(UWORD *)AT.WorkPointer,&num,AN.NoScrat2,&den) ) goto FromNorm;
2378  if ( k < 0 && den < 0 ) {
2379  *AN.NoScrat2 = 1;
2380  den = -1;
2381  if ( AddLong((UWORD *)AT.WorkPointer,num
2382  ,AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2383  goto FromNorm;
2384  }
2385  else if ( k > 0 && den > 0 ) {
2386  *AN.NoScrat2 = 1;
2387  den = 1;
2388  if ( AddLong((UWORD *)AT.WorkPointer,num,
2389  AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2390  goto FromNorm;
2391  }
2392 
2393  }
2394  else if ( *to == -SNUMBER ) { /* No rounding needed */
2395  if ( to[1] < 0 ) { *AT.WorkPointer = -to[1]; num = -1; }
2396  else if ( to[1] == 0 ) goto NormZero;
2397  else { *AT.WorkPointer = to[1]; num = 1; }
2398  }
2399  else goto NoInteg;
2400  if ( num == 0 ) goto NormZero;
2401  ncoef = REDLENG(ncoef);
2402  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,num) )
2403  goto FromNorm;
2404  ncoef = INCLENG(ncoef);
2405  break;
2406  }
2407 NoInteg:;
2408 /*
2409  Fall through if it cannot be resolved
2410 */
2411  default :
2412 defaultcase:;
2413  if ( *t < FUNCTION ) {
2414  MLOCK(ErrorMessageLock);
2415  MesPrint("Illegal code in Norm");
2416 #ifdef DEBUGON
2417  {
2418  UBYTE OutBuf[140];
2419  AO.OutFill = AO.OutputLine = OutBuf;
2420  t = term;
2421  AO.OutSkip = 3;
2422  FiniLine();
2423  i = *t;
2424  while ( --i >= 0 ) {
2425  TalToLine((UWORD)(*t++));
2426  TokenToLine((UBYTE *)" ");
2427  }
2428  AO.OutSkip = 0;
2429  FiniLine();
2430  }
2431 #endif
2432  MUNLOCK(ErrorMessageLock);
2433  goto NormMin;
2434  }
2435  if ( *t == REPLACEMENT ) {
2436  if ( AR.Eside != LHSIDE ) ReplaceVeto--;
2437  pcom[ncom++] = t;
2438  break;
2439  }
2440 /*
2441  if ( *t == AM.termfunnum && t[1] == FUNHEAD+2
2442  && t[FUNHEAD] == -DOLLAREXPRESSION ) termflag++;
2443 */
2444  if ( *t == DUMMYFUN || *t == DUMMYTEN ) {}
2445  else {
2446  if ( *t < (FUNCTION + WILDOFFSET) ) {
2447  if ( ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2448  || ( functions[*t-FUNCTION].minnumargs > 0 ) )
2449  && ( ( t[2] & DIRTYFLAG ) != 0 ) ) {
2450 /*
2451  Number of arguments is bounded. And we have not checked.
2452 */
2453  WORD *ta = t + FUNHEAD, *tb = t + t[1];
2454  int numarg = 0;
2455  while ( ta < tb ) { numarg++; NEXTARG(ta) }
2456  if ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2457  && ( numarg >= functions[*t-FUNCTION].maxnumargs ) )
2458  goto NormZero;
2459  if ( ( functions[*t-FUNCTION].minnumargs > 0 )
2460  && ( numarg < functions[*t-FUNCTION].minnumargs ) )
2461  goto NormZero;
2462  }
2463 doflags:
2464  if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION].tabl == 0 ) ) {
2465  t[2] &= ~DIRTYFLAG;
2466  t[2] |= DIRTYSYMFLAG;
2467  }
2468  if ( functions[*t-FUNCTION].commute ) { pnco[nnco++] = t; }
2469  else { pcom[ncom++] = t; }
2470  }
2471  else {
2472  if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION-WILDOFFSET].tabl == 0 ) ) {
2473  t[2] &= ~DIRTYFLAG;
2474  t[2] |= DIRTYSYMFLAG;
2475  }
2476  if ( functions[*t-FUNCTION-WILDOFFSET].commute ) {
2477  pnco[nnco++] = t;
2478  }
2479  else { pcom[ncom++] = t; }
2480  }
2481  }
2482 
2483  /* Now hunt for contractible indices */
2484 
2485  if ( ( *t < (FUNCTION + WILDOFFSET)
2486  && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) || (
2487  *t >= (FUNCTION + WILDOFFSET)
2488  && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION ) ) {
2489  if ( *t >= GAMMA && *t <= GAMMASEVEN ) t++;
2490  t += FUNHEAD;
2491  while ( t < r ) {
2492  if ( *t == AM.vectorzero ) goto NormZero;
2493  if ( *t >= AM.OffsetIndex && ( *t >= AM.DumInd
2494  || ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2495  pcon[ncon++] = t;
2496  }
2497  else if ( *t == FUNNYWILD ) { t++; }
2498  t++;
2499  }
2500  }
2501  else {
2502  t += FUNHEAD;
2503  while ( t < r ) {
2504  if ( *t > 0 ) {
2505 /*
2506  Here we should worry about a recursion
2507  A problem is the possibility of a construct
2508  like f(mu+nu)
2509 */
2510  t += *t;
2511  }
2512  else if ( *t <= -FUNCTION ) t++;
2513  else if ( *t == -INDEX ) {
2514  if ( t[1] >= AM.OffsetIndex &&
2515  ( t[1] >= AM.DumInd || ( t[1] < AM.WilInd
2516  && indices[t[1]-AM.OffsetIndex].dimension ) ) )
2517  pcon[ncon++] = t+1;
2518  t += 2;
2519  }
2520  else if ( *t == -SYMBOL ) {
2521  if ( t[1] >= MAXPOWER && t[1] < 2*MAXPOWER ) {
2522  *t = -SNUMBER;
2523  t[1] -= MAXPOWER;
2524  }
2525  else if ( t[1] < -MAXPOWER && t[1] > -2*MAXPOWER ) {
2526  *t = -SNUMBER;
2527  t[1] += MAXPOWER;
2528  }
2529  else t += 2;
2530  }
2531  else t += 2;
2532  }
2533  }
2534  break;
2535  }
2536  t = r;
2537 TryAgain:;
2538  } while ( t < m );
2539  if ( ANsc ) {
2540  AN.cTerm = ANsc;
2541  r = t = ANsr; m = ANsm;
2542  ANsc = ANsm = ANsr = 0;
2543  goto conscan;
2544  }
2545 /*
2546  #] First scan :
2547  #[ Easy denominators :
2548 
2549  Easy denominators are denominators that can be replaced by
2550  negative powers of individual subterms. This may add to all
2551  our sublists.
2552 
2553 */
2554  if ( nden ) {
2555  for ( k = 0, i = 0; i < nden; i++ ) {
2556  t = pden[i];
2557  if ( ( t[2] & DIRTYFLAG ) == 0 ) continue;
2558  r = t + t[1]; m = t + FUNHEAD;
2559  if ( m >= r ) {
2560  for ( j = i+1; j < nden; j++ ) pden[j-1] = pden[j];
2561  nden--;
2562  for ( j = 0; j < nnco; j++ ) if ( pnco[j] == t ) break;
2563  for ( j++; j < nnco; j++ ) pnco[j-1] = pnco[j];
2564  nnco--;
2565  i--;
2566  }
2567  else {
2568  NEXTARG(m);
2569  if ( m >= r ) continue;
2570 /*
2571  We have more than one argument. Split the function.
2572 */
2573  if ( k == 0 ) {
2574  k = 1; to = termout; from = term;
2575  }
2576  while ( from < t ) *to++ = *from++;
2577  m = t + FUNHEAD;
2578  while ( m < r ) {
2579  stop = to;
2580  *to++ = DENOMINATOR;
2581  for ( j = 1; j < FUNHEAD; j++ ) *to++ = 0;
2582  if ( *m < -FUNCTION ) *to++ = *m++;
2583  else if ( *m < 0 ) { *to++ = *m++; *to++ = *m++; }
2584  else {
2585  j = *m; while ( --j >= 0 ) *to++ = *m++;
2586  }
2587  stop[1] = WORDDIF(to,stop);
2588  }
2589  from = r;
2590  if ( i == nden - 1 ) {
2591  stop = term + *term;
2592  while ( from < stop ) *to++ = *from++;
2593  i = *termout = WORDDIF(to,termout);
2594  to = term; from = termout;
2595  while ( --i >= 0 ) *to++ = *from++;
2596  goto Restart;
2597  }
2598  }
2599  }
2600  for ( i = 0; i < nden; i++ ) {
2601  t = pden[i];
2602  if ( ( t[2] & DIRTYFLAG ) == 0 ) continue;
2603  t[2] = 0;
2604  if ( t[FUNHEAD] == -SYMBOL ) {
2605  WORD change;
2606  t += FUNHEAD+1;
2607  change = ExtraSymbol(*t,-1,nsym,ppsym,&ncoef);
2608  nsym += change;
2609  ppsym += change * 2;
2610  goto DropDen;
2611  }
2612  else if ( t[FUNHEAD] == -SNUMBER ) {
2613  t += FUNHEAD+1;
2614  if ( *t == 0 ) goto NormInf;
2615  if ( *t < 0 ) { *AT.WorkPointer = -*t; j = -1; }
2616  else { *AT.WorkPointer = *t; j = 1; }
2617  ncoef = REDLENG(ncoef);
2618  if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,j) )
2619  goto FromNorm;
2620  ncoef = INCLENG(ncoef);
2621  goto DropDen;
2622  }
2623  else if ( t[FUNHEAD] == ARGHEAD ) goto NormInf;
2624  else if ( t[FUNHEAD] > 0 && t[FUNHEAD+ARGHEAD] ==
2625  t[FUNHEAD]-ARGHEAD ) {
2626  /* Only one term */
2627  r = t + t[1] - 1;
2628  t += FUNHEAD + ARGHEAD + 1;
2629  j = *r;
2630  m = r - ABS(*r) + 1;
2631  if ( j != 3 || ( ( *m != 1 ) || ( m[1] != 1 ) ) ) {
2632  ncoef = REDLENG(ncoef);
2633  if ( DivRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,REDLENG(j),(UWORD *)n_coef,&ncoef) ) goto FromNorm;
2634  ncoef = INCLENG(ncoef);
2635  j = ABS(j) - 3;
2636  t[-FUNHEAD-ARGHEAD] -= j;
2637  t[-ARGHEAD-1] -= j;
2638  t[-1] -= j;
2639  m[0] = m[1] = 1;
2640  m[2] = 3;
2641  }
2642  while ( t < m ) {
2643  r = t + t[1];
2644  if ( *t == SYMBOL || *t == DOTPRODUCT ) {
2645  k = t[1];
2646  pden[i][1] -= k;
2647  pden[i][FUNHEAD] -= k;
2648  pden[i][FUNHEAD+ARGHEAD] -= k;
2649  m -= k;
2650  stop = m + 3;
2651  tt = to = t;
2652  from = r;
2653  if ( *t == SYMBOL ) {
2654  t += 2;
2655  while ( t < r ) {
2656  WORD change;
2657  change = ExtraSymbol(*t,-t[1],nsym,ppsym,&ncoef);
2658  nsym += change;
2659  ppsym += change * 2;
2660  t += 2;
2661  }
2662  }
2663  else {
2664  t += 2;
2665  while ( t < r ) {
2666  *ppdot++ = *t++;
2667  *ppdot++ = *t++;
2668  *ppdot++ = -*t++;
2669  ndot++;
2670  }
2671  }
2672  while ( to < stop ) *to++ = *from++;
2673  r = tt;
2674  }
2675  t = r;
2676  }
2677  if ( pden[i][1] == 4+FUNHEAD+ARGHEAD ) {
2678 DropDen:
2679  for ( j = 0; j < nnco; j++ ) {
2680  if ( pden[i] == pnco[j] ) {
2681  --nnco;
2682  while ( j < nnco ) {
2683  pnco[j] = pnco[j+1];
2684  j++;
2685  }
2686  break;
2687  }
2688  }
2689  pden[i--] = pden[--nden];
2690  }
2691  }
2692  }
2693  }
2694 /*
2695  #] Easy denominators :
2696  #[ Index Contractions :
2697 */
2698  if ( ndel ) {
2699  t = pdel;
2700  for ( i = 0; i < ndel; i += 2 ) {
2701  if ( t[0] == t[1] ) {
2702  if ( t[0] == EMPTYINDEX ) {}
2703  else if ( *t < AM.OffsetIndex ) {
2704  k = AC.FixIndices[*t];
2705  if ( k < 0 ) { j = -1; k = -k; }
2706  else if ( k > 0 ) j = 1;
2707  else goto NormZero;
2708  goto WithFix;
2709  }
2710  else if ( *t >= AM.DumInd ) {
2711  k = AC.lDefDim;
2712  if ( k ) goto docontract;
2713  }
2714  else if ( *t >= AM.WilInd ) {
2715  k = indices[*t-AM.OffsetIndex-WILDOFFSET].dimension;
2716  if ( k ) goto docontract;
2717  }
2718  else if ( ( k = indices[*t-AM.OffsetIndex].dimension ) != 0 ) {
2719 docontract:
2720  if ( k > 0 ) {
2721  j = 1;
2722 WithFix: shortnum = k;
2723  ncoef = REDLENG(ncoef);
2724  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),j) )
2725  goto FromNorm;
2726  ncoef = INCLENG(ncoef);
2727  }
2728  else {
2729  WORD change;
2730  change = ExtraSymbol((WORD)(-k),(WORD)1,nsym,ppsym,&ncoef);
2731  nsym += change;
2732  ppsym += change * 2;
2733  }
2734  t[1] = pdel[ndel-1];
2735  t[0] = pdel[ndel-2];
2736 HaveCon:
2737  ndel -= 2;
2738  i -= 2;
2739  }
2740  }
2741  else {
2742  if ( *t < AM.OffsetIndex && t[1] < AM.OffsetIndex ) goto NormZero;
2743  j = *t - AM.OffsetIndex;
2744  if ( j >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2745  || ( *t < AM.WilInd && indices[j].dimension ) ) ) {
2746  for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2747  if ( *t == *m ) {
2748  *t = m[1];
2749  *m++ = pdel[ndel-2];
2750  *m = pdel[ndel-1];
2751  goto HaveCon;
2752  }
2753  else if ( *t == m[1] ) {
2754  *t = *m;
2755  *m++ = pdel[ndel-2];
2756  *m = pdel[ndel-1];
2757  goto HaveCon;
2758  }
2759  }
2760  }
2761  j = t[1]-AM.OffsetIndex;
2762  if ( j >= 0 && ( ( t[1] >= AM.DumInd && AC.lDefDim )
2763  || ( t[1] < AM.WilInd && indices[j].dimension ) ) ) {
2764  for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2765  if ( t[1] == *m ) {
2766  t[1] = m[1];
2767  *m++ = pdel[ndel-2];
2768  *m = pdel[ndel-1];
2769  goto HaveCon;
2770  }
2771  else if ( t[1] == m[1] ) {
2772  t[1] = *m;
2773  *m++ = pdel[ndel-2];
2774  *m = pdel[ndel-1];
2775  goto HaveCon;
2776  }
2777  }
2778  }
2779  t += 2;
2780  }
2781  }
2782  if ( ndel > 0 ) {
2783  if ( nvec ) {
2784  t = pdel;
2785  for ( i = 0; i < ndel; i++ ) {
2786  if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2787  ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2788  r = pvec + 1;
2789  for ( j = 1; j < nvec; j += 2 ) {
2790  if ( *r == *t ) {
2791  if ( i & 1 ) {
2792  *r = t[-1];
2793  *t-- = pdel[--ndel];
2794  i -= 2;
2795  }
2796  else {
2797  *r = t[1];
2798  t[1] = pdel[--ndel];
2799  i--;
2800  }
2801  *t-- = pdel[--ndel];
2802  break;
2803  }
2804  r += 2;
2805  }
2806  }
2807  t++;
2808  }
2809  }
2810  if ( ndel > 0 && ncon ) {
2811  t = pdel;
2812  for ( i = 0; i < ndel; i++ ) {
2813  if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2814  ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2815  for ( j = 0; j < ncon; j++ ) {
2816  if ( *pcon[j] == *t ) {
2817  if ( i & 1 ) {
2818  *pcon[j] = t[-1];
2819  *t-- = pdel[--ndel];
2820  i -= 2;
2821  }
2822  else {
2823  *pcon[j] = t[1];
2824  t[1] = pdel[--ndel];
2825  i--;
2826  }
2827  *t-- = pdel[--ndel];
2828  didcontr++;
2829  r = pcon[j];
2830  for ( j = 0; j < nnco; j++ ) {
2831  m = pnco[j];
2832  if ( r > m && r < m+m[1] ) {
2833  m[2] |= DIRTYSYMFLAG;
2834  break;
2835  }
2836  }
2837  for ( j = 0; j < ncom; j++ ) {
2838  m = pcom[j];
2839  if ( r > m && r < m+m[1] ) {
2840  m[2] |= DIRTYSYMFLAG;
2841  break;
2842  }
2843  }
2844  for ( j = 0; j < neps; j++ ) {
2845  m = peps[j];
2846  if ( r > m && r < m+m[1] ) {
2847  m[2] |= DIRTYSYMFLAG;
2848  break;
2849  }
2850  }
2851  break;
2852  }
2853  }
2854  }
2855  t++;
2856  }
2857  }
2858  }
2859  }
2860  if ( nvec ) {
2861  t = pvec + 1;
2862  for ( i = 3; i < nvec; i += 2 ) {
2863  k = *t - AM.OffsetIndex;
2864  if ( k >= 0 && ( ( *t > AM.DumInd && AC.lDefDim )
2865  || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2866  r = t + 2;
2867  for ( j = i; j < nvec; j += 2 ) {
2868  if ( *r == *t ) { /* Another dotproduct */
2869  *ppdot++ = t[-1];
2870  *ppdot++ = r[-1];
2871  *ppdot++ = 1;
2872  ndot++;
2873  *r-- = pvec[--nvec];
2874  *r = pvec[--nvec];
2875  *t-- = pvec[--nvec];
2876  *t-- = pvec[--nvec];
2877  i -= 2;
2878  break;
2879  }
2880  r += 2;
2881  }
2882  }
2883  t += 2;
2884  }
2885  if ( nvec > 0 && ncon ) {
2886  t = pvec + 1;
2887  for ( i = 1; i < nvec; i += 2 ) {
2888  k = *t - AM.OffsetIndex;
2889  if ( k >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2890  || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2891  for ( j = 0; j < ncon; j++ ) {
2892  if ( *pcon[j] == *t ) {
2893  *pcon[j] = t[-1];
2894  *t-- = pvec[--nvec];
2895  *t-- = pvec[--nvec];
2896  r = pcon[j];
2897  pcon[j] = pcon[--ncon];
2898  i -= 2;
2899  for ( j = 0; j < nnco; j++ ) {
2900  m = pnco[j];
2901  if ( r > m && r < m+m[1] ) {
2902  m[2] |= DIRTYSYMFLAG;
2903  break;
2904  }
2905  }
2906  for ( j = 0; j < ncom; j++ ) {
2907  m = pcom[j];
2908  if ( r > m && r < m+m[1] ) {
2909  m[2] |= DIRTYSYMFLAG;
2910  break;
2911  }
2912  }
2913  for ( j = 0; j < neps; j++ ) {
2914  m = peps[j];
2915  if ( r > m && r < m+m[1] ) {
2916  m[2] |= DIRTYSYMFLAG;
2917  break;
2918  }
2919  }
2920  break;
2921  }
2922  }
2923  }
2924  t += 2;
2925  }
2926  }
2927  }
2928 /*
2929  #] Index Contractions :
2930  #[ NonCommuting Functions :
2931 */
2932  m = fillsetexp;
2933  if ( nnco ) {
2934  for ( i = 0; i < nnco; i++ ) {
2935  t = pnco[i];
2936  if ( ( *t >= (FUNCTION+WILDOFFSET)
2937  && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
2938  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2939  && functions[*t-FUNCTION].spec == 0 ) ) {
2940  DoRevert(t,m);
2941  if ( didcontr ) {
2942  r = t + FUNHEAD;
2943  t += t[1];
2944  while ( r < t ) {
2945  if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
2946  *r = -SNUMBER;
2947  didcontr--;
2948  pnco[i][2] |= DIRTYSYMFLAG;
2949  }
2950  NEXTARG(r)
2951  }
2952  }
2953  }
2954  }
2955 
2956  /* First should come the code for function properties. */
2957 
2958  /* First we test for symmetric properties and the DIRTYSYMFLAG */
2959 
2960  for ( i = 0; i < nnco; i++ ) {
2961  t = pnco[i];
2962  if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) && *t != DOLLAREXPRESSION ) {
2963  l = 0; /* to make the compiler happy */
2964  if ( ( *t >= (FUNCTION+WILDOFFSET)
2965  && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
2966  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
2967  && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
2968  if ( *t >= (FUNCTION+WILDOFFSET) ) {
2969  *t -= WILDOFFSET;
2970  j = FullSymmetrize(BHEAD t,l);
2971  *t += WILDOFFSET;
2972  }
2973  else j = FullSymmetrize(BHEAD t,l);
2974  if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
2975  if ( ( j & 2 ) != 0 ) goto NormZero;
2976  if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
2977  }
2978  }
2979  else t[2] &= ~DIRTYSYMFLAG;
2980  }
2981  }
2982 
2983  /* Non commuting functions are then tested for commutation
2984  rules. If needed their order is exchanged. */
2985 
2986  k = nnco - 1;
2987  for ( i = 0; i < k; i++ ) {
2988  j = i;
2989  while ( Commute(pnco[j],pnco[j+1]) ) {
2990  t = pnco[j]; pnco[j] = pnco[j+1]; pnco[j+1] = t;
2991  l = j-1;
2992  while ( l >= 0 && Commute(pnco[l],pnco[l+1]) ) {
2993  t = pnco[l]; pnco[l] = pnco[l+1]; pnco[l+1] = t;
2994  l--;
2995  }
2996  if ( ++j >= k ) break;
2997  }
2998  }
2999 
3000  /* Finally they are written to output. gamma matrices
3001  are bundled if possible */
3002 
3003  for ( i = 0; i < nnco; i++ ) {
3004  t = pnco[i];
3005  if ( *t == IDFUNCTION ) AN.idfunctionflag = 1;
3006  if ( *t >= GAMMA && *t <= GAMMASEVEN ) {
3007  WORD gtype;
3008  to = m;
3009  *m++ = GAMMA;
3010  m++;
3011  FILLFUN(m)
3012  *m++ = stype = t[FUNHEAD]; /* type of string */
3013  j = 0;
3014  nnum = 0;
3015  do {
3016  r = t + t[1];
3017  if ( *t == GAMMAFIVE ) {
3018  gtype = GAMMA5; t += FUNHEAD; goto onegammamatrix; }
3019  else if ( *t == GAMMASIX ) {
3020  gtype = GAMMA6; t += FUNHEAD; goto onegammamatrix; }
3021  else if ( *t == GAMMASEVEN ) {
3022  gtype = GAMMA7; t += FUNHEAD; goto onegammamatrix; }
3023  t += FUNHEAD+1;
3024  while ( t < r ) {
3025  gtype = *t;
3026 onegammamatrix:
3027  if ( gtype == GAMMA5 ) {
3028  if ( j == GAMMA1 ) j = GAMMA5;
3029  else if ( j == GAMMA5 ) j = GAMMA1;
3030  else if ( j == GAMMA7 ) ncoef = -ncoef;
3031  if ( nnum & 1 ) ncoef = -ncoef;
3032  }
3033  else if ( gtype == GAMMA6 || gtype == GAMMA7 ) {
3034  if ( nnum & 1 ) {
3035  if ( gtype == GAMMA6 ) gtype = GAMMA7;
3036  else gtype = GAMMA6;
3037  }
3038  if ( j == GAMMA1 ) j = gtype;
3039  else if ( j == GAMMA5 ) {
3040  j = gtype;
3041  if ( j == GAMMA7 ) ncoef = -ncoef;
3042  }
3043  else if ( j != gtype ) goto NormZero;
3044  else {
3045  shortnum = 2;
3046  ncoef = REDLENG(ncoef);
3047  if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),1) ) goto FromNorm;
3048  ncoef = INCLENG(ncoef);
3049  }
3050  }
3051  else {
3052  *m++ = gtype; nnum++;
3053  }
3054  t++;
3055  }
3056 
3057  } while ( ( ++i < nnco ) && ( *(t = pnco[i]) >= GAMMA
3058  && *t <= GAMMASEVEN ) && ( t[FUNHEAD] == stype ) );
3059  i--;
3060  if ( j ) {
3061  k = WORDDIF(m,to) - FUNHEAD-1;
3062  r = m;
3063  from = m++;
3064  while ( --k >= 0 ) *from-- = *--r;
3065  *from = j;
3066  }
3067  to[1] = WORDDIF(m,to);
3068  }
3069  else if ( *t < 0 ) {
3070  *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3071  FILLFUN3(m)
3072  }
3073  else {
3074  if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3075  && *t != REPLACEMENT && *t != DOLLAREXPRESSION
3076  && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3077  k = t[1];
3078  NCOPY(m,t,k);
3079  }
3080  }
3081 
3082  }
3083 /*
3084  #] NonCommuting Functions :
3085  #[ Commuting Functions :
3086 */
3087  if ( ncom ) {
3088  for ( i = 0; i < ncom; i++ ) {
3089  t = pcom[i];
3090  if ( ( *t >= (FUNCTION+WILDOFFSET)
3091  && functions[*t-FUNCTION-WILDOFFSET].spec == 0 )
3092  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3093  && functions[*t-FUNCTION].spec == 0 ) ) {
3094  DoRevert(t,m);
3095  if ( didcontr ) {
3096  r = t + FUNHEAD;
3097  t += t[1];
3098  while ( r < t ) {
3099  if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3100  *r = -SNUMBER;
3101  didcontr--;
3102  pcom[i][2] |= DIRTYSYMFLAG;
3103  }
3104  NEXTARG(r)
3105  }
3106  }
3107  }
3108  }
3109 
3110  /* Now we test for symmetric properties and the DIRTYSYMFLAG */
3111 
3112  for ( i = 0; i < ncom; i++ ) {
3113  t = pcom[i];
3114  if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) ) {
3115  l = 0; /* to make the compiler happy */
3116  if ( ( *t >= (FUNCTION+WILDOFFSET)
3117  && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3118  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3119  && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3120  if ( *t >= (FUNCTION+WILDOFFSET) ) {
3121  *t -= WILDOFFSET;
3122  j = FullSymmetrize(BHEAD t,l);
3123  *t += WILDOFFSET;
3124  }
3125  else j = FullSymmetrize(BHEAD t,l);
3126  if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3127  if ( ( j & 2 ) != 0 ) goto NormZero;
3128  if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3129  }
3130  }
3131  else t[2] &= ~DIRTYSYMFLAG;
3132  }
3133  }
3134 /*
3135  Sort the functions
3136  From a purists point of view this can be improved.
3137  There arel slow and fast arguments and no conversions are
3138  taken into account here.
3139 */
3140  for ( i = 1; i < ncom; i++ ) {
3141  for ( j = i; j > 0; j-- ) {
3142  WORD jj,kk;
3143  jj = j-1;
3144  t = pcom[jj];
3145  r = pcom[j];
3146  if ( *t < 0 ) {
3147  if ( *r < 0 ) { if ( *t >= *r ) goto NextI; }
3148  else { if ( -*t <= *r ) goto NextI; }
3149  goto jexch;
3150  }
3151  else if ( *r < 0 ) {
3152  if ( *t < -*r ) goto NextI;
3153  goto jexch;
3154  }
3155  else if ( *t != *r ) {
3156  if ( *t < *r ) goto NextI;
3157 jexch: t = pcom[j]; pcom[j] = pcom[jj]; pcom[jj] = t;
3158  continue;
3159  }
3160  if ( AC.properorderflag ) {
3161  if ( ( *t >= (FUNCTION+WILDOFFSET)
3162  && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
3163  || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3164  && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) ) {}
3165  else {
3166  WORD *s1, *s2, *ss1, *ss2;
3167  s1 = t+FUNHEAD; s2 = r+FUNHEAD;
3168  ss1 = t + t[1]; ss2 = r + r[1];
3169  while ( s1 < ss1 && s2 < ss2 ) {
3170  k = CompArg(s1,s2);
3171  if ( k > 0 ) goto jexch;
3172  if ( k < 0 ) goto NextI;
3173  NEXTARG(s1)
3174  NEXTARG(s2)
3175  }
3176  if ( s1 < ss1 ) goto jexch;
3177  goto NextI;
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  else
3192  {
3193  k = t[1] - FUNHEAD;
3194  kk = r[1] - FUNHEAD;
3195  t += FUNHEAD;
3196  r += FUNHEAD;
3197  while ( k > 0 && kk > 0 ) {
3198  if ( *t < *r ) goto NextI;
3199  else if ( *t++ > *r++ ) goto jexch;
3200  k--; kk--;
3201  }
3202  if ( k > 0 ) goto jexch;
3203  goto NextI;
3204  }
3205  }
3206 NextI:;
3207  }
3208  for ( i = 0; i < ncom; i++ ) {
3209  t = pcom[i];
3210  if ( *t == THETA || *t == THETA2 ) {
3211  if ( ( k = DoTheta(BHEAD t) ) == 0 ) goto NormZero;
3212  else if ( k < 0 ) {
3213  k = t[1];
3214  NCOPY(m,t,k);
3215  }
3216  }
3217  else if ( *t == DELTA2 || *t == DELTAP ) {
3218  if ( ( k = DoDelta(t) ) == 0 ) goto NormZero;
3219  else if ( k < 0 ) {
3220  k = t[1];
3221  NCOPY(m,t,k);
3222  }
3223  }
3224  else if ( *t == AR.PolyFunInv && AR.PolyFunType == 2 ) {
3225 /*
3226  If there are two arguments, exchange them, change the
3227  name of the function and go to dealing with PolyRatFun.
3228 */
3229  WORD *mm, *tt = t, numt = 0;
3230  tt += FUNHEAD;
3231  while ( tt < t+t[1] ) { numt++; NEXTARG(tt) }
3232  if ( numt == 2 ) {
3233  tt = t; mm = m; k = t[1];
3234  NCOPY(mm,tt,k)
3235  mm = m+FUNHEAD;
3236  NEXTARG(mm);
3237  tt = t+FUNHEAD;
3238  if ( *mm < 0 ) {
3239  if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3240  else { *tt++ = *mm++; *tt++ = *mm++; }
3241  }
3242  else {
3243  k = *mm; NCOPY(tt,mm,k)
3244  }
3245  mm = m+FUNHEAD;
3246  if ( *mm < 0 ) {
3247  if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3248  else { *tt++ = *mm++; *tt++ = *mm++; }
3249  }
3250  else {
3251  k = *mm; NCOPY(tt,mm,k)
3252  }
3253  *t = AR.PolyFun;
3254  t[2] |= MUSTCLEANPRF;
3255  goto regularratfun;
3256  }
3257  }
3258  else if ( *t == AR.PolyFun ) {
3259  if ( AR.PolyFunType == 1 ) { /* Regular PolyFun with one argument */
3260  if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3261  t[1] == FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) goto NormZero;
3262  if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) {
3263  if ( AN.PolyNormFlag == 0 ) {
3264  AN.PolyNormFlag = 1;
3265  AN.PolyFunTodo = 0;
3266  }
3267  }
3268  k = t[1];
3269  NCOPY(m,t,k);
3270  }
3271  else if ( AR.PolyFunType == 2 ) {
3272 /*
3273  PolyRatFun.
3274  Regular type: Two arguments
3275  Power expanded: One argument. Here to be treated as
3276  AR.PolyFunType == 1, but with power cutoff.
3277 */
3278 regularratfun:;
3279 /*
3280  First check for zeroes.
3281 */
3282  if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3283  t[1] > FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) {
3284  u = t + FUNHEAD + 2;
3285  if ( *u < 0 ) {
3286  if ( *u <= -FUNCTION ) {}
3287  else if ( t[1] == FUNHEAD+4 && t[FUNHEAD+2] == -SNUMBER
3288  && t[FUNHEAD+3] == 0 ) goto NormPRF;
3289  else if ( t[1] == FUNHEAD+4 ) goto NormZero;
3290  }
3291  else if ( t[1] == *u+FUNHEAD+2 ) goto NormZero;
3292  }
3293  else {
3294  u = t+FUNHEAD; NEXTARG(u);
3295  if ( *u == -SNUMBER && u[1] == 0 ) goto NormInf;
3296  }
3297  if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3298  else if ( i < ncom-1 && pcom[i+1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3299  k = t[1];
3300  if ( AN.PolyNormFlag ) {
3301  if ( AR.PolyFunExp == 0 ) {
3302  AN.PolyFunTodo = 0;
3303  NCOPY(m,t,k);
3304  }
3305  else if ( AR.PolyFunExp == 1 ) { /* get highest divergence */
3306  if ( PolyFunMode == 0 ) {
3307  NCOPY(m,t,k);
3308  AN.PolyFunTodo = 1;
3309  }
3310  else {
3311  WORD *mmm = m;
3312  NCOPY(m,t,k);
3313  if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3314  goto FromNorm;
3315  m = mmm+mmm[1];
3316  }
3317  }
3318  else {
3319  if ( PolyFunMode == 0 ) {
3320  NCOPY(m,t,k);
3321  AN.PolyFunTodo = 1;
3322  }
3323  else {
3324  WORD *mmm = m;
3325  NCOPY(m,t,k);
3326  if ( ExpandRat(BHEAD mmm) != 0 )
3327  goto FromNorm;
3328  m = mmm+mmm[1];
3329  }
3330  }
3331  }
3332  else {
3333  if ( AR.PolyFunExp == 0 ) {
3334  AN.PolyFunTodo = 0;
3335  NCOPY(m,t,k);
3336  }
3337  else if ( AR.PolyFunExp == 1 ) { /* get highest divergence */
3338  WORD *mmm = m;
3339  NCOPY(m,t,k);
3340  if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3341  goto FromNorm;
3342  m = mmm+mmm[1];
3343  }
3344  else {
3345  WORD *mmm = m;
3346  NCOPY(m,t,k);
3347  if ( ExpandRat(BHEAD mmm) != 0 )
3348  goto FromNorm;
3349  m = mmm+mmm[1];
3350  }
3351  }
3352  }
3353  }
3354  else if ( *t > 0 ) {
3355  if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3356  && *t != REPLACEMENT && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3357  k = t[1];
3358  NCOPY(m,t,k);
3359  }
3360  else {
3361  *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3362  FILLFUN3(m)
3363  }
3364  }
3365  }
3366 /*
3367  #] Commuting Functions :
3368  #[ Track Replace_ :
3369 */
3370  if ( ReplaceVeto < 0 ) {
3371 /*
3372  We found one (or more) replace_ functions and all other
3373  functions are 'clean' (no dirty flag).
3374  Now we check whether one of these functions can be used.
3375  Thus far the functions go from fillsetexp to m.
3376  Somewhere in there there are -ReplaceVeto occurrences of REPLACEMENT.
3377  Hunt for the first one that fits the bill.
3378  Note that replace_ is a commuting function.
3379 */
3380  WORD *ma = fillsetexp, *mb, *mc;
3381  while ( ma < m ) {
3382  mb = ma + ma[1];
3383  if ( *ma != REPLACEMENT ) {
3384  ma = mb;
3385  continue;
3386  }
3387  if ( *ma == REPLACEMENT && ReplaceType == -1 ) {
3388  mc = ma;
3389  ReplaceType = 0;
3390  if ( AN.RSsize < 2*ma[1]+SUBEXPSIZE ) {
3391  if ( AN.ReplaceScrat ) M_free(AN.ReplaceScrat,"AN.ReplaceScrat");
3392  AN.RSsize = 2*ma[1]+SUBEXPSIZE+40;
3393  AN.ReplaceScrat = (WORD *)Malloc1((AN.RSsize+1)*sizeof(WORD),"AN.ReplaceScrat");
3394  }
3395  ma += FUNHEAD;
3396  ReplaceSub = AN.ReplaceScrat;
3397  ReplaceSub += SUBEXPSIZE;
3398  while ( ma < mb ) {
3399  if ( *ma > 0 ) goto NoRep;
3400  if ( *ma <= -FUNCTION ) {
3401  *ReplaceSub++ = FUNTOFUN;
3402  *ReplaceSub++ = 4;
3403  *ReplaceSub++ = -*ma++;
3404  if ( *ma > -FUNCTION ) goto NoRep;
3405  *ReplaceSub++ = -*ma++;
3406  }
3407  else if ( ma+4 > mb ) goto NoRep;
3408  else {
3409  if ( *ma == -SYMBOL ) {
3410  if ( ma[2] == -SYMBOL && ma+4 <= mb )
3411  *ReplaceSub++ = SYMTOSYM;
3412  else if ( ma[2] == -SNUMBER && ma+4 <= mb ) {
3413  *ReplaceSub++ = SYMTONUM;
3414  if ( ReplaceType == 0 ) {
3415  oldtoprhs = C->numrhs;
3416  oldcpointer = C->Pointer - C->Buffer;
3417  }
3418  ReplaceType = 1;
3419  }
3420  else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3421  *ReplaceSub++ = SYMTONUM;
3422  *ReplaceSub++ = 4;
3423  *ReplaceSub++ = ma[1];
3424  *ReplaceSub++ = 0;
3425  ma += 2+ARGHEAD;
3426  continue;
3427  }
3428 /*
3429  Next is the subexpression. We have to test that
3430  it isn't vector-like or index-like
3431 */
3432  else if ( ma[2] > 0 ) {
3433  WORD *sstop, *ttstop, n;
3434  ss = ma+2;
3435  sstop = ss + *ss;
3436  ss += ARGHEAD;
3437  while ( ss < sstop ) {
3438  tt = ss + *ss;
3439  ttstop = tt - ABS(tt[-1]);
3440  ss++;
3441  while ( ss < ttstop ) {
3442  if ( *ss == INDEX ) goto NoRep;
3443  ss += ss[1];
3444  }
3445  ss = tt;
3446  }
3447  subtype = SYMTOSUB;
3448  if ( ReplaceType == 0 ) {
3449  oldtoprhs = C->numrhs;
3450  oldcpointer = C->Pointer - C->Buffer;
3451  }
3452  ReplaceType = 1;
3453  ss = AddRHS(AT.ebufnum,1);
3454  tt = ma+2;
3455  n = *tt - ARGHEAD;
3456  tt += ARGHEAD;
3457  while ( (ss + n + 10) > C->Top ) ss = DoubleCbuffer(AT.ebufnum,ss,14);
3458  while ( --n >= 0 ) *ss++ = *tt++;
3459  *ss++ = 0;
3460  C->rhs[C->numrhs+1] = ss;
3461  C->Pointer = ss;
3462  *ReplaceSub++ = subtype;
3463  *ReplaceSub++ = 4;
3464  *ReplaceSub++ = ma[1];
3465  *ReplaceSub++ = C->numrhs;
3466  ma += 2 + ma[2];
3467  continue;
3468  }
3469  else goto NoRep;
3470  }
3471  else if ( ( *ma == -VECTOR || *ma == -MINVECTOR ) && ma+4 <= mb ) {
3472  if ( ma[2] == -VECTOR ) {
3473  if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOVEC;
3474  else *ReplaceSub++ = VECTOMIN;
3475  }
3476  else if ( ma[2] == -MINVECTOR ) {
3477  if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOMIN;
3478  else *ReplaceSub++ = VECTOVEC;
3479  }
3480 /*
3481  Next is a vector-like subexpression
3482  Search for vector nature first
3483 */
3484  else if ( ma[2] > 0 ) {
3485  WORD *sstop, *ttstop, *w, *mm, n, count;
3486  WORD *v1, *v2 = 0;
3487  if ( *ma == -MINVECTOR ) {
3488  ss = ma+2;
3489  sstop = ss + *ss;
3490  ss += ARGHEAD;
3491  while ( ss < sstop ) {
3492  ss += *ss;
3493  ss[-1] = -ss[-1];
3494  }
3495  *ma = -VECTOR;
3496  }
3497  ss = ma+2;
3498  sstop = ss + *ss;
3499  ss += ARGHEAD;
3500  while ( ss < sstop ) {
3501  tt = ss + *ss;
3502  ttstop = tt - ABS(tt[-1]);
3503  ss++;
3504  count = 0;
3505  while ( ss < ttstop ) {
3506  if ( *ss == INDEX ) {
3507  n = ss[1] - 2; ss += 2;
3508  while ( --n >= 0 ) {
3509  if ( *ss < MINSPEC ) count++;
3510  ss++;
3511  }
3512  }
3513  else ss += ss[1];
3514  }
3515  if ( count != 1 ) goto NoRep;
3516  ss = tt;
3517  }
3518  subtype = VECTOSUB;
3519  if ( ReplaceType == 0 ) {
3520  oldtoprhs = C->numrhs;
3521  oldcpointer = C->Pointer - C->Buffer;
3522  }
3523  ReplaceType = 1;
3524  mm = AddRHS(AT.ebufnum,1);
3525  *ReplaceSub++ = subtype;
3526  *ReplaceSub++ = 4;
3527  *ReplaceSub++ = ma[1];
3528  *ReplaceSub++ = C->numrhs;
3529  w = ma+2;
3530  n = *w - ARGHEAD;
3531  w += ARGHEAD;
3532  while ( (mm + n + 10) > C->Top )
3533  mm = DoubleCbuffer(AT.ebufnum,mm,15);
3534  while ( --n >= 0 ) *mm++ = *w++;
3535  *mm++ = 0;
3536  C->rhs[C->numrhs+1] = mm;
3537  C->Pointer = mm;
3538  mm = AddRHS(AT.ebufnum,1);
3539  w = ma+2;
3540  n = *w - ARGHEAD;
3541  w += ARGHEAD;
3542  while ( (mm + n + 13) > C->Top )
3543  mm = DoubleCbuffer(AT.ebufnum,mm,16);
3544  sstop = w + n;
3545  while ( w < sstop ) {
3546  tt = w + *w; ttstop = tt - ABS(tt[-1]);
3547  ss = mm; mm++; w++;
3548  while ( w < ttstop ) { /* Subterms */
3549  if ( *w != INDEX ) {
3550  n = w[1];
3551  NCOPY(mm,w,n);
3552  }
3553  else {
3554  v1 = mm;
3555  *mm++ = *w++;
3556  *mm++ = n = *w++;
3557  n -= 2;
3558  while ( --n >= 0 ) {
3559  if ( *w >= MINSPEC ) *mm++ = *w++;
3560  else v2 = w++;
3561  }
3562  n = WORDDIF(mm,v1);
3563  if ( n != v1[1] ) {
3564  if ( n <= 2 ) mm -= 2;
3565  else v1[1] = n;
3566  *mm++ = VECTOR;
3567  *mm++ = 4;
3568  *mm++ = *v2;
3569  *mm++ = FUNNYVEC;
3570  }
3571  }
3572  }
3573  while ( w < tt ) *mm++ = *w++;
3574  *ss = WORDDIF(mm,ss);
3575  }
3576  *mm++ = 0;
3577  C->rhs[C->numrhs+1] = mm;
3578  C->Pointer = mm;
3579  if ( mm > C->Top ) {
3580  MLOCK(ErrorMessageLock);
3581  MesPrint("Internal error in Normalize with extra compiler buffer");
3582  MUNLOCK(ErrorMessageLock);
3583  Terminate(-1);
3584  }
3585  ma += 2 + ma[2];
3586  continue;
3587  }
3588  else goto NoRep;
3589  }
3590  else if ( *ma == -INDEX ) {
3591  if ( ( ma[2] == -INDEX || ma[2] == -VECTOR )
3592  && ma+4 <= mb )
3593  *ReplaceSub++ = INDTOIND;
3594  else if ( ma[1] >= AM.OffsetIndex ) {
3595  if ( ma[2] == -SNUMBER && ma+4 <= mb
3596  && ma[3] >= 0 && ma[3] < AM.OffsetIndex )
3597  *ReplaceSub++ = INDTOIND;
3598  else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3599  *ReplaceSub++ = INDTOIND;
3600  *ReplaceSub++ = 4;
3601  *ReplaceSub++ = ma[1];
3602  *ReplaceSub++ = 0;
3603  ma += 2+ARGHEAD;
3604  continue;
3605  }
3606  else goto NoRep;
3607  }
3608  else goto NoRep;
3609  }
3610  else goto NoRep;
3611  *ReplaceSub++ = 4;
3612  *ReplaceSub++ = ma[1];
3613  *ReplaceSub++ = ma[3];
3614  ma += 4;
3615  }
3616 
3617  }
3618  AN.ReplaceScrat[1] = ReplaceSub-AN.ReplaceScrat;
3619 /*
3620  Success. This means that we have to remove the replace_
3621  from the functions. It starts at mc and end at mb.
3622 */
3623  while ( mb < m ) *mc++ = *mb++;
3624  m = mc;
3625  break;
3626 NoRep:
3627  if ( ReplaceType > 0 ) {
3628  C->numrhs = oldtoprhs;
3629  C->Pointer = C->Buffer + oldcpointer;
3630  }
3631  ReplaceType = -1;
3632  if ( ++ReplaceVeto >= 0 ) break;
3633  }
3634  ma = mb;
3635  }
3636  }
3637 /*
3638  #] Track Replace_ :
3639  #[ LeviCivita tensors :
3640 */
3641  if ( neps ) {
3642  to = m;
3643  for ( i = 0; i < neps; i++ ) { /* Put the indices in order */
3644  t = peps[i];
3645  if ( ( t[2] & DIRTYSYMFLAG ) != DIRTYSYMFLAG ) continue;
3646  t[2] &= ~DIRTYSYMFLAG;
3647  if ( AR.Eside == LHSIDE || AR.Eside == LHSIDEX ) {
3648  /* Potential problems with FUNNYWILD */
3649 /*
3650  First make sure all FUNNIES are at the end.
3651  Then sort separately
3652 */
3653  r = t + FUNHEAD;
3654  m = tt = t + t[1];
3655  while ( r < m ) {
3656  if ( *r != FUNNYWILD ) { r++; continue; }
3657  k = r[1]; u = r + 2;
3658  while ( u < tt ) {
3659  u[-2] = *u;
3660  if ( *u != FUNNYWILD ) ncoef = -ncoef;
3661  u++;
3662  }
3663  tt[-2] = FUNNYWILD; tt[-1] = k; m -= 2;
3664  }
3665  t += FUNHEAD;
3666  do {
3667  for ( r = t + 1; r < m; r++ ) {
3668  if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3669  else if ( *r == *t ) goto NormZero;
3670  }
3671  t++;
3672  } while ( t < m );
3673  do {
3674  for ( r = t + 2; r < tt; r += 2 ) {
3675  if ( r[1] < t[1] ) {
3676  k = r[1]; r[1] = t[1]; t[1] = k; ncoef = -ncoef; }
3677  else if ( r[1] == t[1] ) goto NormZero;
3678  }
3679  t += 2;
3680  } while ( t < tt );
3681  }
3682  else {
3683  m = t + t[1];
3684  t += FUNHEAD;
3685  do {
3686  for ( r = t + 1; r < m; r++ ) {
3687  if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3688  else if ( *r == *t ) goto NormZero;
3689  }
3690  t++;
3691  } while ( t < m );
3692  }
3693  }
3694 
3695  /* Sort the tensors */
3696 
3697  for ( i = 0; i < (neps-1); i++ ) {
3698  t = peps[i];
3699  for ( j = i+1; j < neps; j++ ) {
3700  r = peps[j];
3701  if ( t[1] > r[1] ) {
3702  peps[i] = m = r; peps[j] = r = t; t = m;
3703  }
3704  else if ( t[1] == r[1] ) {
3705  k = t[1] - FUNHEAD;
3706  m = t + FUNHEAD;
3707  r += FUNHEAD;
3708  do {
3709  if ( *r < *m ) {
3710  m = peps[j]; peps[j] = t; peps[i] = t = m;
3711  break;
3712  }
3713  else if ( *r++ > *m++ ) break;
3714  } while ( --k > 0 );
3715  }
3716  }
3717  }
3718  m = to;
3719  for ( i = 0; i < neps; i++ ) {
3720  t = peps[i];
3721  k = t[1];
3722  NCOPY(m,t,k);
3723  }
3724  }
3725 /*
3726  #] LeviCivita tensors :
3727  #[ Delta :
3728 */
3729  if ( ndel ) {
3730  r = t = pdel;
3731  for ( i = 0; i < ndel; i += 2, r += 2 ) {
3732  if ( r[1] < r[0] ) { k = *r; *r = r[1]; r[1] = k; }
3733  }
3734  for ( i = 2; i < ndel; i += 2, t += 2 ) {
3735  r = t + 2;
3736  for ( j = i; j < ndel; j += 2 ) {
3737  if ( *r > *t ) { r += 2; }
3738  else if ( *r < *t ) {
3739  k = *r; *r++ = *t; *t++ = k;
3740  k = *r; *r++ = *t; *t-- = k;
3741  }
3742  else {
3743  if ( *++r < t[1] ) {
3744  k = *r; *r = t[1]; t[1] = k;
3745  }
3746  r++;
3747  }
3748  }
3749  }
3750  t = pdel;
3751  *m++ = DELTA;
3752  *m++ = ndel + 2;
3753  i = ndel;
3754  NCOPY(m,t,i);
3755  }
3756 /*
3757  #] Delta :
3758  #[ Loose Vectors/Indices :
3759 */
3760  if ( nind ) {
3761  t = pind;
3762  for ( i = 0; i < nind; i++ ) {
3763  r = t + 1;
3764  for ( j = i+1; j < nind; j++ ) {
3765  if ( *r < *t ) {
3766  k = *r; *r = *t; *t = k;
3767  }
3768  r++;
3769  }
3770  t++;
3771  }
3772  t = pind;
3773  *m++ = INDEX;
3774  *m++ = nind + 2;
3775  i = nind;
3776  NCOPY(m,t,i);
3777  }
3778 /*
3779  #] Loose Vectors/Indices :
3780  #[ Vectors :
3781 */
3782  if ( nvec ) {
3783  t = pvec;
3784  for ( i = 2; i < nvec; i += 2 ) {
3785  r = t + 2;
3786  for ( j = i; j < nvec; j += 2 ) {
3787  if ( *r == *t ) {
3788  if ( *++r < t[1] ) {
3789  k = *r; *r = t[1]; t[1] = k;
3790  }
3791  r++;
3792  }
3793  else if ( *r < *t ) {
3794  k = *r; *r++ = *t; *t++ = k;
3795  k = *r; *r++ = *t; *t-- = k;
3796  }
3797  else { r += 2; }
3798  }
3799  t += 2;
3800  }
3801  t = pvec;
3802  *m++ = VECTOR;
3803  *m++ = nvec + 2;
3804  i = nvec;
3805  NCOPY(m,t,i);
3806  }
3807 /*
3808  #] Vectors :
3809  #[ Dotproducts :
3810 */
3811  if ( ndot ) {
3812  to = m;
3813  m = t = pdot;
3814  i = ndot;
3815  while ( --i >= 0 ) {
3816  if ( *t > t[1] ) { j = *t; *t = t[1]; t[1] = j; }
3817  t += 3;
3818  }
3819  t = m;
3820  ndot *= 3;
3821  m += ndot;
3822  while ( t < (m-3) ) {
3823  r = t + 3;
3824  do {
3825  if ( *r == *t ) {
3826  if ( *++r == *++t ) {
3827  r++;
3828  if ( ( *r < MAXPOWER && t[1] < MAXPOWER )
3829  || ( *r > -MAXPOWER && t[1] > -MAXPOWER ) ) {
3830  t++;
3831  *t += *r;
3832  if ( *t > MAXPOWER || *t < -MAXPOWER ) {
3833  MLOCK(ErrorMessageLock);
3834  MesPrint("Exponent of dotproduct out of range: %d",*t);
3835  MUNLOCK(ErrorMessageLock);
3836  goto NormMin;
3837  }
3838  ndot -= 3;
3839  *r-- = *--m;
3840  *r-- = *--m;
3841  *r = *--m;
3842  if ( !*t ) {
3843  ndot -= 3;
3844  *t-- = *--m;
3845  *t-- = *--m;
3846  *t = *--m;
3847  t -= 3;
3848  break;
3849  }
3850  }
3851  else if ( *r < *++t ) {
3852  k = *r; *r++ = *t; *t = k;
3853  }
3854  else r++;
3855  t -= 2;
3856  }
3857  else if ( *r < *t ) {
3858  k = *r; *r++ = *t; *t++ = k;
3859  k = *r; *r++ = *t; *t = k;
3860  t -= 2;
3861  }
3862  else { r += 2; t--; }
3863  }
3864  else if ( *r < *t ) {
3865  k = *r; *r++ = *t; *t++ = k;
3866  k = *r; *r++ = *t; *t++ = k;
3867  k = *r; *r++ = *t; *t = k;
3868  t -= 2;
3869  }
3870  else { r += 3; }
3871  } while ( r < m );
3872  t += 3;
3873  }
3874  m = to;
3875  t = pdot;
3876  if ( ( i = ndot ) > 0 ) {
3877  *m++ = DOTPRODUCT;
3878  *m++ = i + 2;
3879  NCOPY(m,t,i);
3880  }
3881  }
3882 /*
3883  #] Dotproducts :
3884  #[ Symbols :
3885 */
3886  if ( nsym ) {
3887  nsym <<= 1;
3888  t = psym;
3889  *m++ = SYMBOL;
3890  r = m;
3891  *m++ = ( i = nsym ) + 2;
3892  if ( i ) { do {
3893  if ( !*t ) {
3894  if ( t[1] < (2*MAXPOWER) ) { /* powers of i */
3895  if ( t[1] & 1 ) { *m++ = 0; *m++ = 1; }
3896  else *r -= 2;
3897  if ( *++t & 2 ) ncoef = -ncoef;
3898  t++;
3899  }
3900  }
3901  else if ( *t <= NumSymbols && *t > -2*MAXPOWER ) { /* Put powers in range */
3902  if ( ( ( ( t[1] > symbols[*t].maxpower ) && ( symbols[*t].maxpower < MAXPOWER ) ) ||
3903  ( ( t[1] < symbols[*t].minpower ) && ( symbols[*t].minpower > -MAXPOWER ) ) ) &&
3904  ( t[1] < 2*MAXPOWER ) && ( t[1] > -2*MAXPOWER ) ) {
3905  if ( i <= 2 || t[2] != *t ) goto NormZero;
3906  }
3907  if ( AN.ncmod == 1 && ( AC.modmode & ALSOPOWERS ) != 0 ) {
3908  if ( AC.cmod[0] == 1 ) t[1] = 0;
3909  else if ( t[1] >= 0 ) t[1] = 1 + (t[1]-1)%(AC.cmod[0]-1);
3910  else {
3911  t[1] = -1 - (-t[1]-1)%(AC.cmod[0]-1);
3912  if ( t[1] < 0 ) t[1] += (AC.cmod[0]-1);
3913  }
3914  }
3915  if ( ( t[1] < (2*MAXPOWER) && t[1] >= MAXPOWER )
3916  || ( t[1] > -(2*MAXPOWER) && t[1] <= -MAXPOWER ) ) {
3917  MLOCK(ErrorMessageLock);
3918  MesPrint("Exponent out of range: %d",t[1]);
3919  MUNLOCK(ErrorMessageLock);
3920  goto NormMin;
3921  }
3922  if ( AT.TrimPower && AR.PolyFunVar == *t && t[1] > AR.PolyFunPow ) {
3923  goto NormZero;
3924  }
3925  else if ( t[1] ) {
3926  *m++ = *t++;
3927  *m++ = *t++;
3928  }
3929  else { *r -= 2; t += 2; }
3930  }
3931  else {
3932  *m++ = *t++; *m++ = *t++;
3933  }
3934  } while ( (i-=2) > 0 ); }
3935  if ( *r <= 2 ) m = r-1;
3936  }
3937 /*
3938  #] Symbols :
3939  #[ Do Replace_ :
3940 */
3941  stop = (WORD *)(((UBYTE *)(termout)) + AM.MaxTer);
3942  i = ABS(ncoef);
3943  if ( ( m + i ) > stop ) {
3944  MLOCK(ErrorMessageLock);
3945  MesPrint("Term too complex during normalization");
3946  MUNLOCK(ErrorMessageLock);
3947  goto NormMin;
3948  }
3949  if ( ReplaceType >= 0 ) {
3950  t = n_coef;
3951  i--;
3952  NCOPY(m,t,i);
3953  *m++ = ncoef;
3954  t = termout;
3955  *t = WORDDIF(m,t);
3956  if ( ReplaceType == 0 ) {
3957  AT.WorkPointer = termout+*termout;
3958  WildFill(BHEAD term,termout,AN.ReplaceScrat);
3959  termout = term + *term;
3960  }
3961  else {
3962  AT.WorkPointer = r = termout + *termout;
3963  WildFill(BHEAD r,termout,AN.ReplaceScrat);
3964  i = *r; m = term;
3965  NCOPY(m,r,i);
3966  termout = m;
3967 
3968 
3969  r = m = term;
3970  r += *term; r -= ABS(r[-1]);
3971  m++;
3972  while ( m < r ) {
3973  if ( *m >= FUNCTION && m[1] > FUNHEAD &&
3974  functions[*m-FUNCTION].spec != TENSORFUNCTION )
3975  m[2] |= DIRTYFLAG;
3976  m += m[1];
3977  }
3978  }
3979 /*
3980  The next 'reset' cannot be done. We still need the expression
3981  in the buffer. Note though that this may cause a runaway pointer
3982  if we are not very careful.
3983 
3984  C->numrhs = oldtoprhs;
3985  C->Pointer = C->Buffer + oldcpointer;
3986 */
3987  TermFree(n_llnum,"n_llnum");
3988  TermFree(n_coef,"NormCoef");
3989  return(1);
3990  }
3991  else {
3992  t = termout;
3993  k = WORDDIF(m,t);
3994  *t = k + i;
3995  m = term;
3996  NCOPY(m,t,k);
3997  i--;
3998  t = n_coef;
3999  NCOPY(m,t,i);
4000  *m++ = ncoef;
4001  }
4002 /*
4003  #] Do Replace_ :
4004  #[ Errors and Finish :
4005 */
4006 RegEnd:
4007  if ( termout < term + *term && termout >= term ) AT.WorkPointer = term + *term;
4008  else AT.WorkPointer = termout;
4009 /*
4010  if ( termflag ) { We have to assign the term to $variable(s)
4011  TermAssign(term);
4012  }
4013 */
4014  TermFree(n_llnum,"n_llnum");
4015  TermFree(n_coef,"NormCoef");
4016  return(regval);
4017 
4018 NormInf:
4019  MLOCK(ErrorMessageLock);
4020  MesPrint("Division by zero during normalization");
4021  MUNLOCK(ErrorMessageLock);
4022  Terminate(-1);
4023 
4024 NormZZ:
4025  MLOCK(ErrorMessageLock);
4026  MesPrint("0^0 during normalization of term");
4027  MUNLOCK(ErrorMessageLock);
4028  Terminate(-1);
4029 
4030 NormPRF:
4031  MLOCK(ErrorMessageLock);
4032  MesPrint("0/0 in polyratfun during normalization of term");
4033  MUNLOCK(ErrorMessageLock);
4034  Terminate(-1);
4035 
4036 NormZero:
4037  *term = 0;
4038  AT.WorkPointer = termout;
4039  TermFree(n_llnum,"n_llnum");
4040  TermFree(n_coef,"NormCoef");
4041  return(regval);
4042 
4043 NormMin:
4044  TermFree(n_llnum,"n_llnum");
4045  TermFree(n_coef,"NormCoef");
4046  return(-1);
4047 
4048 FromNorm:
4049  MLOCK(ErrorMessageLock);
4050  MesCall("Norm");
4051  MUNLOCK(ErrorMessageLock);
4052  TermFree(n_llnum,"n_llnum");
4053  TermFree(n_coef,"NormCoef");
4054  return(-1);
4055 
4056 /*
4057  #] Errors and Finish :
4058 */
4059 }
4060 
4061 /*
4062  #] Normalize :
4063  #[ ExtraSymbol :
4064 */
4065 
4066 WORD ExtraSymbol(WORD sym, WORD pow, WORD nsym, WORD *ppsym, WORD *ncoef)
4067 {
4068  WORD *m, i;
4069  i = nsym;
4070  m = ppsym - 2;
4071  while ( i > 0 ) {
4072  if ( sym == *m ) {
4073  m++;
4074  if ( pow > 2*MAXPOWER || pow < -2*MAXPOWER
4075  || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
4076  MLOCK(ErrorMessageLock);
4077  MesPrint("Illegal wildcard power combination.");
4078  MUNLOCK(ErrorMessageLock);
4079  Terminate(-1);
4080  }
4081  *m += pow;
4082 
4083  if ( ( sym <= NumSymbols && sym > -MAXPOWER )
4084  && ( symbols[sym].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
4085  *m %= symbols[sym].maxpower;
4086  if ( *m < 0 ) *m += symbols[sym].maxpower;
4087  if ( ( symbols[sym].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
4088  if ( ( ( symbols[sym].maxpower & 1 ) == 0 ) &&
4089  ( *m >= symbols[sym].maxpower/2 ) ) {
4090  *m -= symbols[sym].maxpower/2; *ncoef = -*ncoef;
4091  }
4092  }
4093  }
4094 
4095  if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
4096  MLOCK(ErrorMessageLock);
4097  MesPrint("Power overflow during normalization");
4098  MUNLOCK(ErrorMessageLock);
4099  return(-1);
4100  }
4101  if ( !*m ) {
4102  m--;
4103  while ( i < nsym )
4104  { *m = m[2]; m++; *m = m[2]; m++; i++; }
4105  return(-1);
4106  }
4107  return(0);
4108  }
4109  else if ( sym < *m ) {
4110  m -= 2;
4111  i--;
4112  }
4113  else break;
4114  }
4115  m = ppsym;
4116  while ( i < nsym )
4117  { m--; m[2] = *m; m--; m[2] = *m; i++; }
4118  *m++ = sym;
4119  *m = pow;
4120  return(1);
4121 }
4122 
4123 /*
4124  #] ExtraSymbol :
4125  #[ DoTheta :
4126 */
4127 
4128 WORD DoTheta(PHEAD WORD *t)
4129 {
4130  GETBIDENTITY
4131  WORD k, *r1, *r2, *tstop, type;
4132  WORD ia, *ta, *tb, *stopa, *stopb;
4133  if ( AC.BracketNormalize ) return(-1);
4134  type = *t;
4135  k = t[1];
4136  tstop = t + k;
4137  t += FUNHEAD;
4138  if ( k <= FUNHEAD ) return(1);
4139  r1 = t;
4140  NEXTARG(r1)
4141  if ( r1 == tstop ) {
4142 /*
4143  One argument
4144 */
4145  if ( *t == ARGHEAD ) {
4146  if ( type == THETA ) return(1);
4147  else return(0); /* THETA2 */
4148  }
4149  if ( *t < 0 ) {
4150  if ( *t == -SNUMBER ) {
4151  if ( t[1] < 0 ) return(0);
4152  else {
4153  if ( type == THETA2 && t[1] == 0 ) return(0);
4154  else return(1);
4155  }
4156  }
4157  return(-1);
4158  }
4159  k = t[*t-1];
4160  if ( *t == ABS(k)+1+ARGHEAD ) {
4161  if ( k > 0 ) return(1);
4162  else return(0);
4163  }
4164  return(-1);
4165  }
4166 /*
4167  At least two arguments
4168 */
4169  r2 = r1;
4170  NEXTARG(r2)
4171  if ( r2 < tstop ) return(-1); /* More than 2 arguments */
4172 /*
4173  Note now that zero has to be treated specially
4174  We take the criteria from the symmetrize routine
4175 */
4176  if ( *t == -SNUMBER && *r1 == -SNUMBER ) {
4177  if ( t[1] > r1[1] ) return(0);
4178  else if ( t[1] < r1[1] ) {
4179  return(1);
4180  }
4181  else if ( type == THETA ) return(1);
4182  else return(0); /* THETA2 */
4183  }
4184  else if ( t[1] == 0 && *t == -SNUMBER ) {
4185  if ( *r1 > 0 ) { }
4186  else if ( *t < *r1 ) return(1);
4187  else if ( *t > *r1 ) return(0);
4188  }
4189  else if ( r1[1] == 0 && *r1 == -SNUMBER ) {
4190  if ( *t > 0 ) { }
4191  else if ( *t < *r1 ) return(1);
4192  else if ( *t > *r1 ) return(0);
4193  }
4194  r2 = AT.WorkPointer;
4195  if ( *t < 0 ) {
4196  ta = r2;
4197  ToGeneral(t,ta,0);
4198  r2 += *r2;
4199  }
4200  else ta = t;
4201  if ( *r1 < 0 ) {
4202  tb = r2;
4203  ToGeneral(r1,tb,0);
4204  }
4205  else tb = r1;
4206  stopa = ta + *ta;
4207  stopb = tb + *tb;
4208  ta += ARGHEAD; tb += ARGHEAD;
4209  while ( ta < stopa ) {
4210  if ( tb >= stopb ) return(0);
4211  if ( ( ia = CompareTerms(ta,tb,(WORD)1) ) < 0 ) return(0);
4212  if ( ia > 0 ) return(1);
4213  ta += *ta;
4214  tb += *tb;
4215  }
4216  if ( type == THETA ) return(1);
4217  else return(0); /* THETA2 */
4218 }
4219 
4220 /*
4221  #] DoTheta :
4222  #[ DoDelta :
4223 */
4224 
4225 WORD DoDelta(WORD *t)
4226 {
4227  WORD k, *r1, *r2, *tstop, isnum, isnum2, type = *t;
4228  if ( AC.BracketNormalize ) return(-1);
4229  k = t[1];
4230  if ( k <= FUNHEAD ) goto argzero;
4231  if ( k == FUNHEAD+ARGHEAD && t[FUNHEAD] == ARGHEAD ) goto argzero;
4232  tstop = t + k;
4233  t += FUNHEAD;
4234  r1 = t;
4235  NEXTARG(r1)
4236  if ( *t < 0 ) {
4237  k = 1;
4238  if ( *t == -SNUMBER ) { isnum = 1; k = t[1]; }
4239  else isnum = 0;
4240  }
4241  else {
4242  k = t[*t-1];
4243  k = ABS(k);
4244  if ( k == *t-ARGHEAD-1 ) isnum = 1;
4245  else isnum = 0;
4246  k = 1;
4247  }
4248  if ( r1 >= tstop ) { /* Single argument */
4249  if ( !isnum ) return(-1);
4250  if ( k == 0 ) goto argzero;
4251  goto argnonzero;
4252  }
4253  r2 = r1;
4254  NEXTARG(r2)
4255  if ( r2 < tstop ) return(-1);
4256  if ( *r1 < 0 ) {
4257  if ( *r1 == -SNUMBER ) { isnum2 = 1; }
4258  else isnum2 = 0;
4259  }
4260  else {
4261  k = r1[*r1-1];
4262  k = ABS(k);
4263  if ( k == *r1-ARGHEAD-1 ) isnum2 = 1;
4264  else isnum2 = 0;
4265  }
4266  if ( isnum != isnum2 ) return(-1);
4267  tstop = r1;
4268  while ( t < tstop && r1 < r2 ) {
4269  if ( *t != *r1 ) {
4270  if ( !isnum ) return(-1);
4271  goto argnonzero;
4272  }
4273  t++; r1++;
4274  }
4275  if ( t != tstop || r1 != r2 ) {
4276  if ( !isnum ) return(-1);
4277  goto argnonzero;
4278  }
4279 argzero:
4280  if ( type == DELTA2 ) return(1);
4281  else return(0);
4282 argnonzero:
4283  if ( type == DELTA2 ) return(0);
4284  else return(1);
4285 }
4286 
4287 /*
4288  #] DoDelta :
4289  #[ DoRevert :
4290 */
4291 
4292 void DoRevert(WORD *fun, WORD *tmp)
4293 {
4294  WORD *t, *r, *m, *to, *tt, *mm, i, j;
4295  to = fun + fun[1];
4296  r = fun + FUNHEAD;
4297  while ( r < to ) {
4298  if ( *r <= 0 ) {
4299  if ( *r == -REVERSEFUNCTION ) {
4300  m = r; mm = m+1;
4301  while ( mm < to ) *m++ = *mm++;
4302  to--;
4303  (fun[1])--;
4304  fun[2] |= DIRTYSYMFLAG;
4305  }
4306  else if ( *r <= -FUNCTION ) r++;
4307  else {
4308  if ( *r == -INDEX && r[1] < MINSPEC ) *r = -VECTOR;
4309  r += 2;
4310  }
4311  }
4312  else {
4313  if ( ( *r > ARGHEAD )
4314  && ( r[ARGHEAD+1] == REVERSEFUNCTION )
4315  && ( *r == (r[ARGHEAD]+ARGHEAD) )
4316  && ( r[ARGHEAD] == (r[ARGHEAD+2]+4) )
4317  && ( *(r+*r-3) == 1 )
4318  && ( *(r+*r-2) == 1 )
4319  && ( *(r+*r-1) == 3 ) ) {
4320  mm = r;
4321  r += ARGHEAD + 1;
4322  tt = r + r[1];
4323  r += FUNHEAD;
4324  m = tmp;
4325  t = r;
4326  j = 0;
4327  while ( t < tt ) {
4328  NEXTARG(t)
4329  j++;
4330  }
4331  while ( --j >= 0 ) {
4332  i = j;
4333  t = r;
4334  while ( --i >= 0 ) {
4335  NEXTARG(t)
4336  }
4337  if ( *t > 0 ) {
4338  i = *t;
4339  NCOPY(m,t,i);
4340  }
4341  else if ( *t <= -FUNCTION ) *m++ = *t++;
4342  else { *m++ = *t++; *m++ = *t++; }
4343  }
4344  i = WORDDIF(m,tmp);
4345  m = tmp;
4346  t = mm;
4347  r = t + *t;
4348  NCOPY(t,m,i);
4349  m = r;
4350  r = t;
4351  i = WORDDIF(to,m);
4352  NCOPY(t,m,i);
4353  fun[1] = WORDDIF(t,fun);
4354  to = t;
4355  fun[2] |= DIRTYSYMFLAG;
4356  }
4357  else r += *r;
4358  }
4359  }
4360 }
4361 
4362 /*
4363  #] DoRevert :
4364  #] Normalize :
4365  #[ DetCommu :
4366 
4367  Determines the number of terms in an expression that contain
4368  noncommuting objects. This can be used to see whether products of
4369  this expression can be evaluated with binomial coefficients.
4370 
4371  We don't try to be fancy. If a term contains noncommuting objects
4372  we are not looking whether they can commute with complete other
4373  terms.
4374 
4375  If the number gets too large we cut it off.
4376 */
4377 
4378 #define MAXNUMBEROFNONCOMTERMS 2
4379 
4380 WORD DetCommu(WORD *terms)
4381 {
4382  WORD *t, *tnext, *tstop;
4383  WORD num = 0;
4384  if ( *terms == 0 ) return(0);
4385  if ( terms[*terms] == 0 ) return(0);
4386  t = terms;
4387  while ( *t ) {
4388  tnext = t + *t;
4389  tstop = tnext - ABS(tnext[-1]);
4390  t++;
4391  while ( t < tstop ) {
4392  if ( *t >= FUNCTION ) {
4393  if ( functions[*t-FUNCTION].commute ) {
4394  num++;
4395  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4396  break;
4397  }
4398  }
4399  else if ( *t == SUBEXPRESSION ) {
4400  if ( cbuf[t[4]].CanCommu[t[2]] ) {
4401  num++;
4402  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4403  break;
4404  }
4405  }
4406  else if ( *t == EXPRESSION ) {
4407  num++;
4408  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4409  break;
4410  }
4411  else if ( *t == DOLLAREXPRESSION ) {
4412 /*
4413  Technically this is not correct. We have to test first
4414  whether this is MODLOCAL (in TFORM) and if so, use the
4415  local version. Anyway, this should be rare to never
4416  occurring because dollars should be replaced.
4417 */
4418  if ( cbuf[AM.dbufnum].CanCommu[t[2]] ) {
4419  num++;
4420  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4421  break;
4422  }
4423  }
4424  t += t[1];
4425  }
4426  t = tnext;
4427  }
4428  return(num);
4429 }
4430 
4431 /*
4432  #] DetCommu :
4433  #[ DoesCommu :
4434 
4435  Determines the number of noncommuting objects in a term.
4436  If the number gets too large we cut it off.
4437 */
4438 
4439 WORD DoesCommu(WORD *term)
4440 {
4441  WORD *tstop;
4442  WORD num = 0;
4443  if ( *term == 0 ) return(0);
4444  tstop = term + *term;
4445  tstop = tstop - ABS(tstop[-1]);
4446  term++;
4447  while ( term < tstop ) {
4448  if ( ( *term >= FUNCTION ) && ( functions[*term-FUNCTION].commute ) ) {
4449  num++;
4450  if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4451  }
4452  term += term[1];
4453  }
4454  return(num);
4455 }
4456 
4457 /*
4458  #] DoesCommu :
4459  #[ PolyNormPoly :
4460 
4461  Normalizes a polynomial
4462 */
4463 
4464 #ifdef EVALUATEGCD
4465 WORD *PolyNormPoly (PHEAD WORD *Poly) {
4466 
4467  GETBIDENTITY;
4468  WORD *buffer = AT.WorkPointer;
4469  WORD *p;
4470  if ( NewSort(BHEAD0) ) { Terminate(-1); }
4471  AR.CompareRoutine = &CompareSymbols;
4472  while ( *Poly ) {
4473  p = Poly + *Poly;
4474  if ( SymbolNormalize(Poly) < 0 ) return(0);
4475  if ( StoreTerm(BHEAD Poly) ) {
4476  AR.CompareRoutine = &Compare1;
4477  LowerSortLevel();
4478  Terminate(-1);
4479  }
4480  Poly = p;
4481  }
4482  if ( EndSort(BHEAD buffer,1) < 0 ) {
4483  AR.CompareRoutine = &Compare1;
4484  Terminate(-1);
4485  }
4486  p = buffer;
4487  while ( *p ) p += *p;
4488  AR.CompareRoutine = &Compare1;
4489  AT.WorkPointer = p + 1;
4490  return(buffer);
4491 }
4492 #endif
4493 
4494 /*
4495  #] PolyNormPoly :
4496  #[ EvaluateGcd :
4497 
4498  Try to evaluate the GCDFUNCTION gcd_.
4499  This function can have a number of arguments which can be numbers
4500  and/or polynomials. If there are objects that aren't SYMBOLS or numbers
4501  it cannot work currently.
4502 
4503  To make this work properly we have to intervene in proces.c
4504  proces.c: if ( Normalize(BHEAD m) ) {
4505 1060 proces.c: if ( Normalize(BHEAD r) ) {
4506 1126?proces.c: if ( Normalize(BHEAD term) ) {
4507  proces.c: if ( Normalize(BHEAD AT.WorkPointer) ) goto PasErr;
4508 2308!proces.c: if ( ( retnorm = Normalize(BHEAD term) ) != 0 ) {
4509  proces.c: ReNumber(BHEAD term); Normalize(BHEAD term);
4510  proces.c: if ( Normalize(BHEAD v) ) Terminate(-1);
4511  proces.c: if ( Normalize(BHEAD w) ) { LowerSortLevel(); goto PolyCall; }
4512  proces.c: if ( Normalize(BHEAD term) ) goto PolyCall;
4513 */
4514 #ifdef EVALUATEGCD
4515 
4516 WORD *EvaluateGcd(PHEAD WORD *subterm)
4517 {
4518  GETBIDENTITY
4519  WORD *oldworkpointer = AT.WorkPointer, *work1, *work2, *work3;
4520  WORD *t, *tt, *ttt, *t1, *t2, *t3, *t4, *tstop;
4521  WORD ct, nnum;
4522  UWORD gcdnum, stor;
4523  WORD *lnum=n_llnum+1;
4524  WORD *num1, *num2, *num3, *den1, *den2, *den3;
4525  WORD sizenum1, sizenum2, sizenum3, sizeden1, sizeden2, sizeden3;
4526  int i, isnumeric = 0, numarg = 0 /*, sizearg */;
4527  LONG size;
4528 /*
4529  Step 1: Look for -SNUMBER or -SYMBOL arguments.
4530  If encountered, treat everybody with it.
4531 */
4532  tt = subterm + subterm[1]; t = subterm + FUNHEAD;
4533 
4534  while ( t < tt ) {
4535  numarg++;
4536  if ( *t == -SNUMBER ) {
4537  if ( t[1] == 0 ) {
4538 gcdzero:;
4539  MLOCK(ErrorMessageLock);
4540  MesPrint("Trying to take the GCD involving a zero term.");
4541  MUNLOCK(ErrorMessageLock);
4542  return(0);
4543  }
4544  gcdnum = ABS(t[1]);
4545  t1 = subterm + FUNHEAD;
4546  while ( gcdnum > 1 && t1 < tt ) {
4547  if ( *t1 == -SNUMBER ) {
4548  stor = ABS(t1[1]);
4549  if ( stor == 0 ) goto gcdzero;
4550  if ( GcdLong(BHEAD (UWORD *)&stor,1,(UWORD *)&gcdnum,1,
4551  (UWORD *)lnum,&nnum) ) goto FromGCD;
4552  gcdnum = lnum[0];
4553  t1 += 2;
4554  continue;
4555  }
4556  else if ( *t1 == -SYMBOL ) goto gcdisone;
4557  else if ( *t1 < 0 ) goto gcdillegal;
4558 /*
4559  Now we have to go through all the terms in the argument.
4560  This includes long numbers.
4561 */
4562  ttt = t1 + *t1;
4563  ct = *ttt; *ttt = 0;
4564  if ( t1[1] != 0 ) { /* First normalize the argument */
4565  t1 = PolyNormPoly(BHEAD t1+ARGHEAD);
4566  }
4567  else t1 += ARGHEAD;
4568  while ( *t1 ) {
4569  t1 += *t1;
4570  i = ABS(t1[-1]);
4571  t2 = t1 - i;
4572  i = (i-1)/2;
4573  t3 = t2+i-1;
4574  while ( t3 > t2 && *t3 == 0 ) { t3--; i--; }
4575  if ( GcdLong(BHEAD (UWORD *)t2,(WORD)i,(UWORD *)&gcdnum,1,
4576  (UWORD *)lnum,&nnum) ) {
4577  *ttt = ct;
4578  goto FromGCD;
4579  }
4580  gcdnum = lnum[0];
4581  if ( gcdnum == 1 ) {
4582  *ttt = ct;
4583  goto gcdisone;
4584  }
4585  }
4586  *ttt = ct;
4587  t1 = ttt;
4588  AT.WorkPointer = oldworkpointer;
4589  }
4590  if ( gcdnum == 1 ) goto gcdisone;
4591  oldworkpointer[0] = 4;
4592  oldworkpointer[1] = gcdnum;
4593  oldworkpointer[2] = 1;
4594  oldworkpointer[3] = 3;
4595  oldworkpointer[4] = 0;
4596  AT.WorkPointer = oldworkpointer + 5;
4597  return(oldworkpointer);
4598  }
4599  else if ( *t == -SYMBOL ) {
4600  t1 = subterm + FUNHEAD;
4601  i = t[1];
4602  while ( t1 < tt ) {
4603  if ( *t1 == -SNUMBER ) goto gcdisone;
4604  if ( *t1 == -SYMBOL ) {
4605  if ( t1[1] != i ) goto gcdisone;
4606  t1 += 2; continue;
4607  }
4608  if ( *t1 < 0 ) goto gcdillegal;
4609  ttt = t1 + *t1;
4610  ct = *ttt; *ttt = 0;
4611  if ( t1[1] != 0 ) { /* First normalize the argument */
4612  t2 = PolyNormPoly(BHEAD t1+ARGHEAD);
4613  }
4614  else t2 = t1 + ARGHEAD;
4615  while ( *t2 ) {
4616  t3 = t2+1;
4617  t2 = t2 + *t2;
4618  tstop = t2 - ABS(t2[-1]);
4619  while ( t3 < tstop ) {
4620  if ( *t3 != SYMBOL ) {
4621  *ttt = ct;
4622  goto gcdillegal;
4623  }
4624  t4 = t3 + 2;
4625  t3 += t3[1];
4626  while ( t4 < t3 ) {
4627  if ( *t4 == i && t4[1] > 0 ) goto nextterminarg;
4628  t4 += 2;
4629  }
4630  }
4631  *ttt = ct;
4632  goto gcdisone;
4633 nextterminarg:;
4634  }
4635  *ttt = ct;
4636  t1 = ttt;
4637  AT.WorkPointer = oldworkpointer;
4638  }
4639  oldworkpointer[0] = 8;
4640  oldworkpointer[1] = SYMBOL;
4641  oldworkpointer[2] = 4;
4642  oldworkpointer[3] = t[1];
4643  oldworkpointer[4] = 1;
4644  oldworkpointer[5] = 1;
4645  oldworkpointer[6] = 1;
4646  oldworkpointer[7] = 3;
4647  oldworkpointer[8] = 0;
4648  AT.WorkPointer = oldworkpointer+9;
4649  return(oldworkpointer);
4650  }
4651  else if ( *t < 0 ) {
4652 gcdillegal:;
4653  MLOCK(ErrorMessageLock);
4654  MesPrint("Illegal object in gcd_ function. Object not a number or a symbol.");
4655  MUNLOCK(ErrorMessageLock);
4656  goto FromGCD;
4657  }
4658  else if ( ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4659  else if ( t[1] != 0 ) {
4660  ttt = t + *t; ct = *ttt; *ttt = 0;
4661  t = PolyNormPoly(BHEAD t+ARGHEAD);
4662  *ttt = ct;
4663  if ( t[*t] == 0 && ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4664  AT.WorkPointer = oldworkpointer;
4665  t = ttt;
4666  }
4667  t += *t;
4668  }
4669 /*
4670  At this point there are only generic arguments.
4671  There are however still two cases:
4672  1: There is an argument that is purely numerical
4673  In that case we have to take the gcd of all coefficients
4674  2: All arguments are nontrivial polynomials.
4675  Here we don't worry so much about the factor. (???)
4676  We know whether case 1 occurs when isnumeric > 0.
4677  We can look up numarg to get a good starting value.
4678 */
4679  AT.WorkPointer = oldworkpointer;
4680  if ( isnumeric ) {
4681  t = subterm + FUNHEAD;
4682  for ( i = 1; i < isnumeric; i++ ) {
4683  NEXTARG(t);
4684  }
4685  if ( t[1] != 0 ) { /* First normalize the argument */
4686  ttt = t + *t; ct = *ttt; *ttt = 0;
4687  t = PolyNormPoly(BHEAD t+ARGHEAD);
4688  *ttt = ct;
4689  }
4690  t += *t;
4691  i = (ABS(t[-1])-1)/2;
4692  den1 = t - 1 - i;
4693  num1 = den1 - i;
4694  sizenum1 = sizeden1 = i;
4695  while ( sizenum1 > 1 && num1[sizenum1-1] == 0 ) sizenum1--;
4696  while ( sizeden1 > 1 && den1[sizeden1-1] == 0 ) sizeden1--;
4697  work1 = AT.WorkPointer+1; work2 = work1+sizenum1;
4698  for ( i = 0; i < sizenum1; i++ ) work1[i] = num1[i];
4699  for ( i = 0; i < sizeden1; i++ ) work2[i] = den1[i];
4700  num1 = work1; den1 = work2;
4701  AT.WorkPointer = work2 = work2 + sizeden1;
4702  t = subterm + FUNHEAD;
4703  while ( t < tt ) {
4704  ttt = t + *t; ct = *ttt; *ttt = 0;
4705  if ( t[1] != 0 ) {
4706  t = PolyNormPoly(BHEAD t+ARGHEAD);
4707  }
4708  else t += ARGHEAD;
4709  while ( *t ) {
4710  t += *t;
4711  i = (ABS(t[-1])-1)/2;
4712  den2 = t - 1 - i;
4713  num2 = den2 - i;
4714  sizenum2 = sizeden2 = i;
4715  while ( sizenum2 > 1 && num2[sizenum2-1] == 0 ) sizenum2--;
4716  while ( sizeden2 > 1 && den2[sizeden2-1] == 0 ) sizeden2--;
4717  num3 = AT.WorkPointer;
4718  if ( GcdLong(BHEAD (UWORD *)num2,sizenum2,(UWORD *)num1,sizenum1,
4719  (UWORD *)num3,&sizenum3) ) goto FromGCD;
4720  sizenum1 = sizenum3;
4721  for ( i = 0; i < sizenum1; i++ ) num1[i] = num3[i];
4722  den3 = AT.WorkPointer;
4723  if ( GcdLong(BHEAD (UWORD *)den2,sizeden2,(UWORD *)den1,sizeden1,
4724  (UWORD *)den3,&sizeden3) ) goto FromGCD;
4725  sizeden1 = sizeden3;
4726  for ( i = 0; i < sizeden1; i++ ) den1[i] = den3[i];
4727  if ( sizenum1 == 1 && num1[0] == 1 && sizeden1 == 1 && den1[1] == 1 )
4728  goto gcdisone;
4729  }
4730  *ttt = ct;
4731  t = ttt;
4732  AT.WorkPointer = work2;
4733  }
4734  AT.WorkPointer = oldworkpointer;
4735 /*
4736  Now copy the GCD to the 'output'
4737 */
4738  if ( sizenum1 > sizeden1 ) {
4739  while ( sizenum1 > sizeden1 ) den1[sizeden1++] = 0;
4740  }
4741  else if ( sizenum1 < sizeden1 ) {
4742  while ( sizenum1 < sizeden1 ) num1[sizenum1++] = 0;
4743  }
4744  t = oldworkpointer;
4745  i = 2*sizenum1+1;
4746  *t++ = i+1;
4747  if ( num1 != t ) { NCOPY(t,num1,sizenum1); }
4748  else t += sizenum1;
4749  if ( den1 != t ) { NCOPY(t,den1,sizeden1); }
4750  else t += sizeden1;
4751  *t++ = i;
4752  *t++ = 0;
4753  AT.WorkPointer = t;
4754  return(oldworkpointer);
4755  }
4756 /*
4757  Now the real stuff with only polynomials.
4758  Pick up the shortest term to start.
4759  We are a bit brutish about this.
4760 */
4761  t = subterm + FUNHEAD;
4762  AT.WorkPointer += AM.MaxTer/sizeof(WORD);
4763  work2 = AT.WorkPointer;
4764 /*
4765  sizearg = subterm[1];
4766 */
4767  i = 0; work3 = 0;
4768  while ( t < tt ) {
4769  i++;
4770  work1 = AT.WorkPointer;
4771  ttt = t + *t; ct = *ttt; *ttt = 0;
4772  t = PolyNormPoly(BHEAD t+ARGHEAD);
4773  if ( *work1 < AT.WorkPointer-work1 ) {
4774 /*
4775  sizearg = AT.WorkPointer-work1;
4776 */
4777  numarg = i;
4778  work3 = work1;
4779  }
4780  *ttt = ct; t = ttt;
4781  }
4782  *AT.WorkPointer++ = 0;
4783 /*
4784  We have properly normalized arguments and the shortest is indicated in work3
4785 */
4786  work1 = work3;
4787  while ( *work2 ) {
4788  if ( work2 != work3 ) {
4789  work1 = PolyGCD2(BHEAD work1,work2);
4790  }
4791  while ( *work2 ) work2 += *work2;
4792  work2++;
4793  }
4794  work2 = work1;
4795  while ( *work2 ) work2 += *work2;
4796  size = work2 - work1 + 1;
4797  t = oldworkpointer;
4798  NCOPY(t,work1,size);
4799  AT.WorkPointer = t;
4800  return(oldworkpointer);
4801 
4802 gcdisone:;
4803  oldworkpointer[0] = 4;
4804  oldworkpointer[1] = 1;
4805  oldworkpointer[2] = 1;
4806  oldworkpointer[3] = 3;
4807  oldworkpointer[4] = 0;
4808  AT.WorkPointer = oldworkpointer+5;
4809  return(oldworkpointer);
4810 FromGCD:
4811  MLOCK(ErrorMessageLock);
4812  MesCall("EvaluateGcd");
4813  MUNLOCK(ErrorMessageLock);
4814  return(0);
4815 }
4816 
4817 #endif
4818 
4819 /*
4820  #] EvaluateGcd :
4821  #[ TreatPolyRatFun :
4822 
4823  if ( AR.PolyFunExp == 1 ) we have to trim the contents of the polyratfun
4824  down to its most divergent term and give it coefficient +1. This is done
4825  by taking the terms with the least power in the variable in the numerator
4826  and in the denominator and then combine them.
4827  Answer is either PolyRatFun(ep^n,1) or PolyRatFun(1,1) or PolyRatFun(1,ep^n)
4828 */
4829 
4830 int TreatPolyRatFun(PHEAD WORD *prf)
4831 {
4832  WORD *t, *tstop, *r, *rstop, *m, *mstop;
4833  WORD exp1 = MAXPOWER, exp2 = MAXPOWER;
4834  t = prf+FUNHEAD;
4835  if ( *t < 0 ) {
4836  if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
4837  if ( exp1 > 1 ) exp1 = 1;
4838  t += 2;
4839  }
4840  else {
4841  if ( exp1 > 0 ) exp1 = 0;
4842  NEXTARG(t)
4843  }
4844  }
4845  else {
4846  tstop = t + *t;
4847  t += ARGHEAD;
4848  while ( t < tstop ) {
4849 /*
4850  Now look for the minimum power of AR.PolyFunVar
4851 */
4852  r = t+1;
4853  t += *t;
4854  rstop = t - ABS(t[-1]);
4855  while ( r < rstop ) {
4856  if ( *r != SYMBOL ) { r += r[1]; continue; }
4857  m = r;
4858  mstop = m + m[1];
4859  m += 2;
4860  while ( m < mstop ) {
4861  if ( *m == AR.PolyFunVar ) {
4862  if ( m[1] < exp1 ) exp1 = m[1];
4863  break;
4864  }
4865  m += 2;
4866  }
4867  if ( m == mstop ) {
4868  if ( exp1 > 0 ) exp1 = 0;
4869  }
4870  break;
4871  }
4872  if ( r == rstop ) {
4873  if ( exp1 > 0 ) exp1 = 0;
4874  }
4875  }
4876  t = tstop;
4877  }
4878  if ( *t < 0 ) {
4879  if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
4880  if ( exp2 > 1 ) exp2 = 1;
4881  }
4882  else {
4883  if ( exp2 > 0 ) exp2 = 0;
4884  }
4885  }
4886  else {
4887  tstop = t + *t;
4888  t += ARGHEAD;
4889  while ( t < tstop ) {
4890 /*
4891  Now look for the minimum power of AR.PolyFunVar
4892 */
4893  r = t+1;
4894  t += *t;
4895  rstop = t - ABS(t[-1]);
4896  while ( r < rstop ) {
4897  if ( *r != SYMBOL ) { r += r[1]; continue; }
4898  m = r;
4899  mstop = m + m[1];
4900  m += 2;
4901  while ( m < mstop ) {
4902  if ( *m == AR.PolyFunVar ) {
4903  if ( m[1] < exp2 ) exp2 = m[1];
4904  break;
4905  }
4906  m += 2;
4907  }
4908  if ( m == mstop ) {
4909  if ( exp2 > 0 ) exp2 = 0;
4910  }
4911  break;
4912  }
4913  if ( r == rstop ) {
4914  if ( exp2 > 0 ) exp2 = 0;
4915  }
4916  }
4917  }
4918 /*
4919  Now we can compose the output.
4920  Notice that the output can never be longer than the input provided
4921  we never can have arguments that consist of just a function.
4922 */
4923  exp1 = exp1-exp2;
4924 /* if ( exp1 > 0 ) exp1 = 0; */
4925  t = prf+FUNHEAD;
4926  if ( exp1 == 0 ) {
4927  *t++ = -SNUMBER; *t++ = 1;
4928  *t++ = -SNUMBER; *t++ = 1;
4929  }
4930  else if ( exp1 > 0 ) {
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  *t++ = -SNUMBER; *t++ = 1;
4942  }
4943  else {
4944  *t++ = -SNUMBER; *t++ = 1;
4945  if ( exp1 == -1 ) {
4946  *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
4947  }
4948  else {
4949  *t++ = 8+ARGHEAD;
4950  *t++ = 0;
4951  FILLARG(t);
4952  *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
4953  *t++ = -exp1; *t++ = 1; *t++ = 1; *t++ = 3;
4954  }
4955  }
4956  prf[2] = 0; /* Clean */
4957  prf[1] = t - prf;
4958  return(0);
4959 }
4960 
4961 /*
4962  #] TreatPolyRatFun :
4963  #[ DropCoefficient :
4964 */
4965 
4966 void DropCoefficient(PHEAD WORD *term)
4967 {
4968  GETBIDENTITY
4969  WORD *t = term + *term;
4970  WORD n, na;
4971  n = t[-1]; na = ABS(n);
4972  t -= na;
4973  if ( n == 3 && t[0] == 1 && t[1] == 1 ) return;
4974  *AN.RepPoint = 1;
4975  t[0] = 1; t[1] = 1; t[2] = 3;
4976  *term -= (na-3);
4977 }
4978 
4979 /*
4980  #] DropCoefficient :
4981  #[ DropSymbols :
4982 */
4983 
4984 void DropSymbols(PHEAD WORD *term)
4985 {
4986  GETBIDENTITY
4987  WORD *tend = term + *term, *t1, *t2, *tstop;
4988  tstop = tend - ABS(tend[-1]);
4989  t1 = term+1;
4990  while ( t1 < tstop ) {
4991  if ( *t1 == SYMBOL ) {
4992  *AN.RepPoint = 1;
4993  t2 = t1+t1[1];
4994  while ( t2 < tend ) *t1++ = *t2++;
4995  *term = t1 - term;
4996  break;
4997  }
4998  t1 += t1[1];
4999  }
5000 }
5001 
5002 /*
5003  #] DropSymbols :
5004  #[ SymbolNormalize :
5005 */
5014 int SymbolNormalize(WORD *term)
5015 {
5016  GETIDENTITY
5017  WORD buffer[7*NORMSIZE], *t, *b, *bb, *tt, *m, *tstop;
5018  int i;
5019  b = buffer;
5020  *b++ = SYMBOL; *b++ = 2;
5021  t = term + *term;
5022  tstop = t - ABS(t[-1]);
5023  t = term + 1;
5024  while ( t < tstop ) { /* Step 1: collect symbols */
5025  if ( *t == SYMBOL && t < tstop ) {
5026  for ( i = 2; i < t[1]; i += 2 ) {
5027  bb = buffer+2;
5028  while ( bb < b ) {
5029  if ( bb[0] == t[i] ) { /* add powers */
5030  bb[1] += t[i+1];
5031  if ( bb[1] > MAXPOWER || bb[1] < -MAXPOWER ) {
5032  MLOCK(ErrorMessageLock);
5033  MesPrint("Power in SymbolNormalize out of range");
5034  MUNLOCK(ErrorMessageLock);
5035  return(-1);
5036  }
5037  if ( bb[1] == 0 ) {
5038  b -= 2;
5039  while ( bb < b ) {
5040  bb[0] = bb[2]; bb[1] = bb[3]; bb += 2;
5041  }
5042  }
5043  goto Nexti;
5044  }
5045  else if ( bb[0] > t[i] ) { /* insert it */
5046  m = b;
5047  while ( m > bb ) { m[1] = m[-1]; m[0] = m[-2]; m -= 2; }
5048  b += 2;
5049  bb[0] = t[i];
5050  bb[1] = t[i+1];
5051  goto Nexti;
5052  }
5053  bb += 2;
5054  }
5055  if ( bb >= b ) { /* add it to the end */
5056  *b++ = t[i]; *b++ = t[i+1];
5057  }
5058 Nexti:;
5059  }
5060  }
5061  else {
5062  MLOCK(ErrorMessageLock);
5063  MesPrint("Illegal term in SymbolNormalize");
5064  MUNLOCK(ErrorMessageLock);
5065  return(-1);
5066  }
5067  t += t[1];
5068  }
5069  buffer[1] = b - buffer;
5070 /*
5071  Veto negative powers
5072 */
5073  if ( AT.LeaveNegative == 0 ) {
5074  b = buffer; bb = b + b[1]; b += 3;
5075  while ( b < bb ) {
5076  if ( *b < 0 ) {
5077  MLOCK(ErrorMessageLock);
5078  MesPrint("Negative power in SymbolNormalize");
5079  MUNLOCK(ErrorMessageLock);
5080  return(-1);
5081  }
5082  b += 2;
5083  }
5084  }
5085 /*
5086  Now we use the fact that the new term will not be longer than the old one
5087  Actually it should be shorter when there is more than one subterm!
5088  Copy back.
5089 */
5090  i = buffer[1];
5091  b = buffer; tt = term + 1;
5092  if ( i > 2 ) { NCOPY(tt,b,i) }
5093  if ( tt < tstop ) {
5094  i = term[*term-1];
5095  if ( i < 0 ) i = -i;
5096  *term -= (tstop-tt);
5097  NCOPY(tt,tstop,i)
5098  }
5099  return(0);
5100 }
5101 
5102 /*
5103  #] SymbolNormalize :
5104  #[ TestFunFlag :
5105 
5106  Tests whether a function still has unsubstituted subexpressions
5107  This function has its dirtyflag on!
5108 */
5109 
5110 int TestFunFlag(PHEAD WORD *tfun)
5111 {
5112  WORD *t, *tstop, *r, *rstop, *m, *mstop;
5113  if ( functions[*tfun-FUNCTION].spec ) return(0);
5114  tstop = tfun + tfun[1];
5115  t = tfun + FUNHEAD;
5116  while ( t < tstop ) {
5117  if ( *t < 0 ) { NEXTARG(t); continue; }
5118  rstop = t + *t;
5119  if ( t[1] == 0 ) { t = rstop; continue; }
5120  r = t + ARGHEAD;
5121  while ( r < rstop ) { /* Here we loop over terms */
5122  m = r+1; mstop = r+*r; mstop -= ABS(mstop[-1]);
5123  while ( m < mstop ) { /* Loop over the subterms */
5124  if ( *m == SUBEXPRESSION || *m == EXPRESSION || *m == DOLLAREXPRESSION ) return(1);
5125  if ( ( *m >= FUNCTION ) && ( ( m[2] & DIRTYFLAG ) == DIRTYFLAG )
5126  && ( *m != REPLACEMENT ) && TestFunFlag(BHEAD m) ) return(1);
5127  m += m[1];
5128  }
5129  r += *r;
5130  }
5131  t += *t;
5132  }
5133  return(0);
5134 }
5135 
5136 /*
5137  #] TestFunFlag :
5138  #[ BracketNormalize :
5139 */
5140 
5141 WORD BracketNormalize(PHEAD WORD *term)
5142 {
5143  WORD *stop = term+*term-3, *t, *tt, *tstart, *r;
5144  WORD *oldwork = AT.WorkPointer;
5145  WORD *termout;
5146  WORD i, ii, j;
5147  termout = AT.WorkPointer = term+*term;
5148 /*
5149  First collect all functions and sort them
5150 */
5151  tt = termout+1; t = term+1;
5152  while ( t < stop ) {
5153  if ( *t >= FUNCTION ) { i = t[1]; NCOPY(tt,t,i); }
5154  else t += t[1];
5155  }
5156  if ( tt > termout+1 && tt-termout-1 > termout[2] ) { /* sorting */
5157  r = termout+1; ii = tt-r;
5158  for ( i = 0; i < ii-FUNHEAD; i += FUNHEAD ) { /* Bubble sort */
5159  for ( j = i+FUNHEAD; j > 0; j -= FUNHEAD ) {
5160  if ( functions[r[j-FUNHEAD]-FUNCTION].commute
5161  && functions[r[j]-FUNCTION].commute == 0 ) break;
5162  if ( r[j-FUNHEAD] > r[j] ) EXCH(r[j-FUNHEAD],r[j])
5163  else break;
5164  }
5165  }
5166  }
5167 
5168  tstart = tt; t = term + 1; *tt++ = DELTA; *tt++ = 2;
5169  while ( t < stop ) {
5170  if ( *t == DELTA ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5171  else t += t[1];
5172  }
5173  if ( tstart[1] > 2 ) {
5174  for ( r = tstart+2; r < tstart+tstart[1]; r += 2 ) {
5175  if ( r[0] > r[1] ) EXCH(r[0],r[1])
5176  }
5177  }
5178  if ( tstart[1] > 4 ) { /* sorting */
5179  r = tstart+2; ii = tstart[1]-2;
5180  for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5181  for ( j = i+2; j > 0; j -= 2 ) {
5182  if ( r[j-2] > r[j] ) {
5183  EXCH(r[j-2],r[j])
5184  EXCH(r[j-1],r[j+1])
5185  }
5186  else if ( r[j-2] < r[j] ) break;
5187  else {
5188  if ( r[j-1] > r[j+1] ) EXCH(r[j-1],r[j+1])
5189  else break;
5190  }
5191  }
5192  }
5193  tt = tstart+tstart[1];
5194  }
5195  else if ( tstart[1] == 2 ) { tt = tstart; }
5196  else tt = tstart+4;
5197 
5198  tstart = tt; t = term + 1; *tt++ = INDEX; *tt++ = 2;
5199  while ( t < stop ) {
5200  if ( *t == INDEX ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5201  else t += t[1];
5202  }
5203  if ( tstart[1] >= 4 ) { /* sorting */
5204  r = tstart+2; ii = tstart[1]-2;
5205  for ( i = 0; i < ii-1; i += 1 ) { /* Bubble sort */
5206  for ( j = i+1; j > 0; j -= 1 ) {
5207  if ( r[j-1] > r[j] ) EXCH(r[j-1],r[j])
5208  else break;
5209  }
5210  }
5211  tt = tstart+tstart[1];
5212  }
5213  else if ( tstart[1] == 2 ) { tt = tstart; }
5214  else tt = tstart+3;
5215 
5216  tstart = tt; t = term + 1; *tt++ = DOTPRODUCT; *tt++ = 2;
5217  while ( t < stop ) {
5218  if ( *t == DOTPRODUCT ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5219  else t += t[1];
5220  }
5221  if ( tstart[1] > 5 ) { /* sorting */
5222  r = tstart+2; ii = tstart[1]-2;
5223  for ( i = 0; i < ii; i += 3 ) {
5224  if ( r[i] > r[i+1] ) EXCH(r[i],r[i+1])
5225  }
5226  for ( i = 0; i < ii-3; i += 3 ) { /* Bubble sort */
5227  for ( j = i+3; j > 0; j -= 3 ) {
5228  if ( r[j-3] < r[j] ) break;
5229  if ( r[j-3] > r[j] ) {
5230  EXCH(r[j-3],r[j])
5231  EXCH(r[j-2],r[j+1])
5232  }
5233  else {
5234  if ( r[j-2] > r[j+1] ) EXCH(r[j-2],r[j+1])
5235  else break;
5236  }
5237  }
5238  }
5239  tt = tstart+tstart[1];
5240  }
5241  else if ( tstart[1] == 2 ) { tt = tstart; }
5242  else {
5243  if ( tstart[2] > tstart[3] ) EXCH(tstart[2],tstart[3])
5244  tt = tstart+5;
5245  }
5246 
5247  tstart = tt; t = term + 1; *tt++ = SYMBOL; *tt++ = 2;
5248  while ( t < stop ) {
5249  if ( *t == SYMBOL ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5250  else t += t[1];
5251  }
5252  if ( tstart[1] > 4 ) { /* sorting */
5253  r = tstart+2; ii = tstart[1]-2;
5254  for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5255  for ( j = i+2; j > 0; j -= 2 ) {
5256  if ( r[j-2] > r[j] ) EXCH(r[j-2],r[j])
5257  else break;
5258  }
5259  }
5260  tt = tstart+tstart[1];
5261  }
5262  else if ( tstart[1] == 2 ) { tt = tstart; }
5263  else tt = tstart+4;
5264 
5265  tstart = tt; t = term + 1; *tt++ = SETSET; *tt++ = 2;
5266  while ( t < stop ) {
5267  if ( *t == SETSET ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5268  else t += t[1];
5269  }
5270  if ( tstart[1] > 4 ) { /* sorting */
5271  r = tstart+2; ii = tstart[1]-2;
5272  for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5273  for ( j = i+2; j > 0; j -= 2 ) {
5274  if ( r[j-2] > r[j] ) {
5275  EXCH(r[j-2],r[j])
5276  EXCH(r[j-1],r[j+1])
5277  }
5278  else break;
5279  }
5280  }
5281  tt = tstart+tstart[1];
5282  }
5283  else if ( tstart[1] == 2 ) { tt = tstart; }
5284  else tt = tstart+4;
5285  *tt++ = 1; *tt++ = 1; *tt++ = 3;
5286  t = term; i = *termout = tt - termout; tt = termout;
5287  NCOPY(t,tt,i);
5288  AT.WorkPointer = oldwork;
5289  return(0);
5290 }
5291 
5292 /*
5293  #] BracketNormalize :
5294 */
WORD Compare1(WORD *, WORD *, WORD)
Definition: sort.c:2536
WORD CompareSymbols(WORD *, WORD *, WORD)
Definition: sort.c:2976
int SymbolNormalize(WORD *term)
Definition: normal.c:5014
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:4333
WORD ** rhs
Definition: structs.h:943
VOID LowerSortLevel()
Definition: sort.c:4727
WORD * Buffer
Definition: structs.h:939
WORD NewSort(PHEAD0)
Definition: sort.c:592
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:682
WORD * AddRHS(int num, int type)
Definition: comtool.c:214