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