FORM  4.3
opera.c
Go to the documentation of this file.
1 
9 /* #[ License : */
10 /*
11  * Copyright (C) 1984-2022 J.A.M. Vermaseren
12  * When using this file you are requested to refer to the publication
13  * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
14  * This is considered a matter of courtesy as the development was paid
15  * for by FOM the Dutch physics granting agency and we would like to
16  * be able to track its scientific use to convince FOM of its value
17  * for the community.
18  *
19  * This file is part of FORM.
20  *
21  * FORM is free software: you can redistribute it and/or modify it under the
22  * terms of the GNU General Public License as published by the Free Software
23  * Foundation, either version 3 of the License, or (at your option) any later
24  * version.
25  *
26  * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
27  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
28  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
29  * details.
30  *
31  * You should have received a copy of the GNU General Public License along
32  * with FORM. If not, see <http://www.gnu.org/licenses/>.
33  */
34 /* #] License : */
35 /*
36  #[ Includes : opera.c
37 */
38 
39 #include "form3.h"
40 /*
41  int hulp;
42 */
43 
44 /*
45  #] Includes :
46  #[ Operations :
47  #[ EpfFind : WORD EpfFind(term,params)
48 
49  Searches for a pair of Levi-Civita tensors that should be
50  contracted.
51  If a match is found its settings are recorded in AT.TMout.
52  type indicates the number of indices that is searched for,
53  unless all are searched for (type = 0).
54  number is the number of tensors that should survive contraction.
55 
56 */
57 
58 WORD EpfFind(PHEAD WORD *term, WORD *params)
59 {
60  GETBIDENTITY
61  WORD *t, *m, *r, n1 = 0, n2, min = -1, count, fac;
62  WORD *c1 = 0, *c2 = 0, sgn = 1;
63  WORD *tstop, *mstop;
64  UWORD *facto = (UWORD *)AT.WorkPointer;
65  WORD ncoef,nfac;
66  WORD number, type;
67  if ( ( AT.WorkPointer = (WORD *)(facto + AM.MaxTal) ) > AT.WorkTop ) {
68  MLOCK(ErrorMessageLock);
69  MesWork();
70  MUNLOCK(ErrorMessageLock);
71  return(-1);
72  }
73  number = params[3];
74  type = params[4];
75  t = term;
76  GETSTOP(t,tstop);
77  t++;
78  if ( !type ) {
79  while ( *t != LEVICIVITA && t < tstop ) t += t[1];
80  if ( t >= tstop ) return(0);
81  m = t;
82  while ( *m == LEVICIVITA && m < tstop ) { n1++; m += m[1]; }
83 AllLev:
84  if ( n1 <= (number+1) || n1 <= 1 ) return(0);
85  mstop = m;
86  m = t + t[1];
87  do {
88  while ( m[1] == t[1] ) {
89  m += FUNHEAD;
90  r = t+FUNHEAD;
91  count = fac = n1 = n2 = t[1] - FUNHEAD;
92  while ( n1 && n2 ) {
93  if ( *m > *r ) {
94  r++; n2--;
95  }
96  else if ( *m < *r ) {
97  m++; n1--;
98  }
99  else {
100  if ( *m >= AM.OffsetIndex &&
101  ( ( *m >= AM.IndDum && AC.lDefDim == fac ) ||
102  ( *m < AM.IndDum &&
103  indices[*m-AM.OffsetIndex].dimension == fac ) ) ) {
104  count--;
105  }
106  r++; m++; n1--; n2--;
107  }
108  }
109  m += n1;
110  if ( min < 0 || count < min ) {
111  c1 = t;
112  c2 = m - fac - FUNHEAD;
113  min = count;
114  }
115  if ( m >= mstop ) break;
116  }
117  t += t[1];
118  } while ( ( m = t + t[1] ) < mstop );
119  }
120  else {
121  fac = type + FUNHEAD;
122  while ( *t != LEVICIVITA && t < tstop ) t += t[1];
123  while ( *t == LEVICIVITA && t < tstop && t[1] != fac ) t += t[1];
124  if ( t >= tstop ) return(0);
125  m = t;
126  while ( *m == LEVICIVITA && m < tstop && m[1] == fac ) { n1++; m += m[1]; }
127  goto AllLev;
128  }
129 /*
130  We have now the two tensors that give the minimum contraction
131  in c1 and c2.
132  Prepare the AT.TMout array;
133 */
134  if ( min < 0 ) return(0); /* No matching pair! */
135  t = c1;
136  mstop = c2;
137  fac = t[1] - FUNHEAD;
138  m = AT.TMout;
139  *m++ = 3 + (min*2); /* The full length */
140  *m++ = CONTRACT;
141  if ( number < 0 ) *m++ = 1;
142  else *m++ = 0;
143  n1 = n2 = t[1] - FUNHEAD;
144  r = c1 + FUNHEAD;
145  c1 = m;
146  m = c2 + FUNHEAD;
147  c2 = c1 + min;
148  while ( n1 && n2 ) {
149  if ( *m > *r ) { *c1++ = *r++; n2--; }
150  else if ( *m < *r ) { *c2++ = *m++; n1--; }
151  else {
152  if ( *m < AM.OffsetIndex || ( *m < AM.IndDum &&
153  ( indices[*m-AM.OffsetIndex].dimension != fac ) ) ||
154  ( *m >= AM.IndDum && AC.lDefDim != fac ) ) {
155  *c1++ = *r++; *c2++ = *m++;
156  }
157  else { if ( ( n1 ^ n2 ) & 1 ) sgn = -sgn; r++; m++; }
158  n1--; n2--;
159  }
160  }
161  if ( n1 ) { NCOPY(c2,m,n1); }
162  else if ( n2 ) { NCOPY(c1,r,n2); }
163  fac -= min;
164  m = t + t[1];
165  while ( m < mstop ) *t++ = *m++;
166  m += m[1];
167  while ( m < tstop ) *t++ = *m++;
168  *t++ = SUBEXPRESSION;
169  *t++ = SUBEXPSIZE;
170  *t++ = -1;
171  *t++ = 1;
172  *t++ = DUMMYBUFFER;
173  FILLSUB(t)
174  r = term;
175  r += *r - 1;
176  mstop = r;
177  ncoef = REDLENG(*r);
178  tstop = t;
179  while ( m < mstop ) *t++ = *m++;
180  if ( Factorial(BHEAD fac,facto,&nfac) || Mully(BHEAD (UWORD *)tstop,&ncoef,facto,nfac) ) {
181  MLOCK(ErrorMessageLock);
182  MesCall("EpfFind");
183  MUNLOCK(ErrorMessageLock);
184  SETERROR(-1)
185  }
186  tstop += (ABS(ncoef))*2;
187  if ( sgn < 0 ) ncoef = -ncoef;
188  ncoef *= 2;
189  *tstop++ = (ncoef<0)?(ncoef-1):(ncoef+1);
190  *term = WORDDIF(tstop,term);
191  return(1);
192 }
193 
194 /*
195  #] EpfFind :
196  #[ EpfCon : WORD EpfCon(term,params,num,level)
197 
198  Contraction of two strings of indices/vectors. They come
199  from LeviCivita tensors that are being contracted.
200  term is the term with the subterm to be replaced.
201  params is the full indicator:
202  Length, number, raise, parameters.
203  Length is the length of the parameter field.
204  number is the number of the operation.
205  raise tells whether level should be raised afterwards.
206  the parameters are the two strings.
207  level is the id level.
208  The factorial has been multiplied in already.
209 
210 */
211 
212 WORD EpfCon(PHEAD WORD *term, WORD *params, WORD num, WORD level)
213 {
214  GETBIDENTITY
215  WORD *kron, *perm, *termout, *tstop, size2;
216  WORD *m, *t, sizes, sgn = 0, i;
217  sizes = *params - 3;
218  kron = AT.WorkPointer;
219  perm = (AT.WorkPointer += sizes);
220  termout = (AT.WorkPointer += sizes);
221  AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
222  if ( AT.WorkPointer > AT.WorkTop ) {
223  MLOCK(ErrorMessageLock);
224  MesWork();
225  MUNLOCK(ErrorMessageLock);
226  return(-1);
227  }
228  params += 2;
229  if ( !(*params++) ) level--;
230  size2 = sizes>>1;
231  if ( !size2 ) goto DoOnce;
232  while ( ( sgn = EpfGen(size2,params,kron,perm,sgn) ) != 0 ) {
233 DoOnce:
234  t = term;
235  GETSTOP(t,tstop);
236  m = termout;
237  tstop -= SUBEXPSIZE;
238  while ( t < tstop ) *m++ = *t++;
239  if ( t[2] != num || *t != SUBEXPRESSION ) {
240  MLOCK(ErrorMessageLock);
241  MesPrint("Serious error in EpfCon");
242  MUNLOCK(ErrorMessageLock);
243  return(-1);
244  }
245  tstop += SUBEXPSIZE;
246  if ( sizes ) {
247  *m++ = DELTA;
248  *m++ = sizes + 2;
249  t = kron;
250  i = sizes;
251  if ( i ) { NCOPY(m,t,i); }
252  }
253  t = tstop;
254  tstop = term + *term;
255  while ( t < tstop ) *m++ = *t++;
256  *termout = WORDDIF(m,termout);
257  m--;
258  if ( sgn < 0 ) *m = - *m;
259  if ( *termout ) {
260  *AN.RepPoint = 1;
261  AR.expchanged = 1;
262  AT.WorkPointer = termout + *termout;
263  if ( Generator(BHEAD termout,level) < 0 ) goto EpfCall;
264  }
265  }
266  AT.WorkPointer = kron;
267  return(0);
268 EpfCall:
269  if ( AM.tracebackflag ) {
270  MLOCK(ErrorMessageLock);
271  MesCall("EpfCon");
272  MUNLOCK(ErrorMessageLock);
273  }
274  return(-1);
275 }
276 
277 /*
278  #] EpfCon :
279  #[ EpfGen : WORD EpfGen(number,inlist,kron,perm,sgn)
280 */
281 
282 WORD EpfGen(WORD number, WORD *inlist, WORD *kron, WORD *perm, WORD sgn)
283 {
284  WORD i, *in2, k, a;
285  if ( !sgn ) {
286  in2 = inlist + number;
287  number *= 2;
288  for ( i = 1; i < number; i += 2 ) {
289  *perm++ = i;
290  *perm++ = i;
291  *kron++ = *inlist++;
292  *kron++ = *in2++;
293  }
294  if ( number <= 0 ) return(0);
295  else return(1);
296  }
297  number *= 2;
298  i = number - 1;
299  while ( ( i -= 2 ) >= 0 ) {
300  if ( ( k = perm[i] ) != i ) {
301  sgn = -sgn;
302  a = kron[i];
303  kron[i] = kron[k];
304  kron[k] = a;
305  }
306  if ( ( k = ( perm[i] += 2 ) ) < number ) {
307  a = kron[i];
308  kron[i] = kron[k];
309  kron[k] = a;
310  sgn = - sgn;
311  for ( k = i + 2; k < number; k += 2 ) perm[k] = k;
312  return(sgn);
313  }
314  }
315  return(0);
316 }
317 
318 /*
319  #] EpfGen :
320  #[ Trick : WORD Trick(in,t)
321 
322  This routine implements the identity:
323  g_(j,mu)*g_(j,nu)*g_(j,ro)=e_(mu,nu,ro,si)*g5_(j)*g_(j,si)
324  +d_(mu,nu)*g_(j,ro)-d_(mu,ro)*g_(j,nu)+d_(nu,ro)*g_(j,mu)
325  which is for 4 dimensions only!
326 
327  Note that z->gamm = 1 if there is no g5 present.
328 
329 */
330 
331 WORD Trick(WORD *in, TRACES *t)
332 {
333  WORD n, n1;
334  n = t->stap;
335  n1 = t->step1;
336  switch ( t->eers[n] ) {
337  case 0: {
338  WORD *p;
339  p = t->pepf + t->mepf[n];
340  *p++ = *in++;
341  *p++ = *in++;
342  *p++ = *in;
343  *p = ++(t->mdum);
344  (t->mepf[n1]) += 4;
345  *in = t->mdum;
346  t->gamm = - t->gamm;
347  t->eers[n] = 5;
348  break;
349  }
350  case 4: {
351  WORD *p;
352  p = t->pdel + t->mdel[n];
353  (t->mepf[n1]) -= 4;
354  (t->mdum)--;
355  *p++ = *in++;
356  *p = *in++;
357  *in = *(t->pepf + t->mepf[n] + 2);
358  (t->mdel[n1]) += 2;
359  t->gamm = - t->gamm;
360  break;
361  }
362  case 3: {
363  t->sign1 = - t->sign1;
364  *(t->pdel + t->mdel[n] + 1) = in[2];
365  in[2] = in[1];
366  break;
367  }
368  case 2: {
369  t->sign1 = - t->sign1;
370  in[2] = in[0];
371  *(t->pdel + t->mdel[n]) = in[1];
372  break;
373  }
374  case 1: {
375  in[2] = *(t->pdel + t->mdel[n] + 1);
376  (t->mdel[n1]) -= 2;
377  break;
378  }
379  default: {
380  return(0);
381  }
382  }
383  return(--(t->eers[n]));
384 }
385 
386 /*
387  #] Trick :
388  #[ Trace4no : WORD Trace4no(number,kron,t)
389 
390  Takes the trace of a string of gamma matrices in 4 dimensions.
391  There is no test for indices or vectors that are the same.
392  The four dimensions refer to the contraction in the algebra:
393  g_(i,a,b,c) = e_(a,b,c,d)*g_(i,5_,d) + g_(i,a)*d_(b,c)
394  - g_(i,b)*d_(a,c) + g_(i,c)*d_(a,b)
395  This simplifies life very much and leads to shorter expressions
396  than in the n dimensional case.
397 
398  Parameters:
399  number: the number of vectors/indices in inlist.
400  inlist: the indices/vectors in the string.
401  kron: the output delta's.
402  gamma5: the potential gamma5 in front.
403  t: the struct for scratch manipulations.
404  stack: the space to put all scratch arrays in.
405 
406  The return value is zero if there are no more terms, 1 if a
407  term was generated with a positive sign and -1 if a term was
408  generated with a negative sign.
409  The value of one is increased to two if the first 4 values
410  in kron should be interpreted as a Levi-Civita tensor.
411 
412  Note that kron should have more places reserved than the number
413  of indices in inlist, because it will contain dummy indices
414  temporarily. In principle there can be 'number*1/4' extra dummies.
415 
416 */
417 
418 WORD Trace4no(WORD number, WORD *kron, TRACES *t)
419 {
420  WORD i;
421  WORD *p, *m;
422  WORD retval, *stop, oldsign;
423  if ( !t->sgn ) { /* Startup */
424  if ( ( number < 0 ) || ( number & 1 ) ) return(0);
425  if ( number <= 2 ) {
426  if ( t->gamma5 == GAMMA5 ) return(0);
427  if ( number == 2 ) {
428  *kron++ = *t->inlist;
429  *kron++ = t->inlist[1];
430  }
431  return(1);
432  }
433  t->sgn = 1;
434  {
435  WORD nhalf = number >> 1;
436  WORD ndouble = number * 2;
437  p = t->eers;
438  t->eers = p; p += nhalf;
439  t->mepf = p; p += nhalf;
440  t->mdel = p; p += nhalf;
441  t->pdel = p; p += number + nhalf;
442  t->pepf = p; p += ndouble;
443  t->e4 = p; p += number;
444  t->e3 = p; p += ndouble;
445  t->nt3 = p; p += nhalf;
446  t->nt4 = p; p += nhalf;
447  t->j3 = p; p += ndouble;
448  t->j4 = p;
449  }
450  t->mepf[0] = 0;
451  t->mdel[0] = 0;
452  t->mdum = AM.mTraceDum;
453  t->kstep = -2;
454  t->step1 = 0;
455  t->sign1 = 1;
456  t->lc3 = -1;
457  t->lc4 = -1;
458  t->gamm = 1;
459 
460  do {
461  t->stap = (t->step1)++;
462  t->kstep += 2;
463  t->eers[t->stap] = 0;
464  t->mepf[t->step1] = t->mepf[t->stap];
465  t->mdel[t->step1] = t->mdel[t->stap];
466 CallTrick: while ( !Trick(t->inlist+t->kstep,t) ) {
467  t->kstep -= 2;
468  t->step1 = (t->stap)--;
469  if ( t->stap < 0 ) {
470  return(0);
471  }
472  }
473  } while ( t->kstep < (number-4) );
474 /*
475  Take now the trace of the leftover matrices.
476  If gamma5 causes the term to vanish there will be a
477  renewed call to Trick for its next term.
478 */
479  t->sign2 = t->sign1;
480  if ( ( t->gamma5 == GAMMA7 ) && ( t->gamm == -1 ) ) {
481  t->sign2 = - t->sign2;
482  }
483  else if ( ( t->gamma5 == GAMMA5 ) && ( t->gamm == 1 ) ) {
484  goto CallTrick;
485  }
486  else if ( ( t->gamma5 == GAMMA1 ) && ( t->gamm == -1 ) ) {
487  goto CallTrick;
488  }
489  p = t->pdel + t->mdel[t->step1];
490  *p++ = t->inlist[t->kstep+2];
491  *p++ = t->inlist[t->kstep+3];
492 /*
493  Now the trace has been expressed in terms of Levi-Civita tensors
494  and Kronecker delta's.
495  The Levi-Civita tensors are in t->pepf
496  and there are t->mepf[step1] elements in this array.
497  The Kronecker delta's are in t->pdel
498  and there are t->mdel[step1] elements in this array.
499 
500  Next we rake the Levi-Civita tensors together such that there
501  is an optimal use of the contractions.
502 */
503  {
504  WORD ae;
505  ae = t->mepf[t->step1];
506  t->ad = t->mdel[t->step1]+2;
507  t->a4 = 0;
508  t->a3 = 0;
509  while ( ( ae -= 4 ) >= 0 ) {
510  if ( t->pepf[ae] > AM.mTraceDum && t->pepf[ae] <= t->mdum ) {
511  p = t->e3 + t->a3;
512  m = t->pepf + ae;
513  for ( i = 0; i < 3; i++ ) {
514  p[3] = m[3-i];
515  *p++ = m[i-4];
516  }
517  t->a3 += 6;
518  ae -= 4;
519  }
520  else {
521  p = t->e4 + t->a4;
522  m = t->pepf + ae;
523  for ( i = 0; i < 4; i++ ) *p++ = *m++;
524  t->a4 += 4;
525  }
526  }
527  }
528 /*
529  Now e3 contains pairs of LeviCivita tensors that have
530  three indices each and a3 is the total number of indices.
531  e4 contains individual tensors with 4 indices.
532  Some indices may be contracted with Kronecker delta's.
533 
534  Contract the e3 tensors first.
535 */
536 
537  while ( t->a3 > 0 ) {
538  t->nt3[++(t->lc3)] = 0;
539  while ( ( t->nt3[t->lc3] = EpfGen(3,t->e3+t->a3-6,
540  t->pdel+t->ad,t->j3+6*t->lc3,oldsign = t->nt3[t->lc3]) ) == 0 ) {
541  if ( oldsign < 0 ) t->sign2 = - t->sign2;
542  (t->lc3)--;
543 NextE3: if ( t->lc3 < 0 ) goto CallTrick;
544  t->ad -= 6;
545  t->a3 += 6;
546  }
547  if ( oldsign ) {
548  if ( oldsign != t->nt3[t->lc3] ) t->sign2 = - t->sign2;
549  }
550  else if ( t->nt3[t->lc3] < 0 ) t->sign2 = - t->sign2;
551  t->a3 -= 6;
552  t->ad += 6;
553  }
554 /*
555  Contract the e4 tensors.
556 */
557  while ( t->a4 > 4 ) {
558  t->nt4[++(t->lc4)] = 0;
559  while ( ( t->nt4[t->lc4] = EpfGen(4,t->e4+t->a4-8,
560  t->pdel+t->ad,t->j4+8*t->lc4,oldsign = t->nt4[t->lc4]) ) == 0 ) {
561  if ( oldsign < 0 ) t->sign2 = - t->sign2;
562  (t->lc4)--;
563 NextE4: if ( t->lc4 < 0 ) goto NextE3;
564  t->ad -= 8;
565  t->a4 += 8;
566  }
567  if ( oldsign ) {
568  if ( oldsign != t->nt4[t->lc4] ) t->sign2 = - t->sign2;
569  }
570  else if ( t->nt4[t->lc4] < 0 ) t->sign2 = - t->sign2;
571  t->a4 -= 8;
572  t->ad += 8;
573  }
574 /*
575  Finally the extra dummy indices can be eliminated.
576  Note that there are currently t->ad words in t->pdel forming
577  t->ad / 2 Kronecker delta's. We are however not allowed to
578  alter anything in these arrays, so the results should be
579  copied to kron.
580 */
581  m = kron;
582  if ( t->a4 == 4 ) {
583  p = t->e4;
584  *m++ = *p++; *m++ = *p++; *m++ = *p++; *m++ = *p++;
585  retval = 2;
586  }
587  else retval = 1;
588  if ( t->sign2 < 0 ) retval = - retval;
589  p = t->pdel;
590  for ( i = 0; i < t->ad; i++ ) *m++ = *p++;
591  p = kron;
592  if ( t->a4 == 4 ) {
593 /*
594  Test for dummies in the last position of the e_.
595 */
596  stop = p + t->ad + 4;
597  p += 3;
598  while ( *p >= AM.mTraceDum && *p <= t->mdum ) {
599  m = p + 1;
600  do {
601  if ( *m == *p ) {
602  *p = m[1];
603  *m = *--stop;
604  m[1] = *--stop;
605  break;
606  }
607  else if ( m[1] == *p ) {
608  *p = *m;
609  *m = *--stop;
610  m[1] = *--stop;
611  break;
612  }
613  else m += 2;
614  } while ( m < stop );
615  }
616  p++;
617  }
618  else stop = p + t->ad;
619  while ( p < (stop-2) ) {
620  while ( *p >= AM.mTraceDum && *p <= t->mdum ) {
621  m = p + 2;
622  do {
623  if ( *m == *p ) {
624  *p = m[1];
625  *m = *--stop;
626  m[1] = *--stop;
627  break;
628  }
629  else if ( m[1] == *p ) {
630  *p = *m;
631  *m = *--stop;
632  m[1] = *--stop;
633  break;
634  }
635  else m += 2;
636  } while ( m < stop );
637  }
638  p++;
639  while ( *p >= AM.mTraceDum && *p <= t->mdum ) {
640  m = p + 1;
641  do {
642  if ( *m == *p ) {
643  *p = m[1];
644  *m = *--stop;
645  m[1] = *--stop;
646  break;
647  }
648  else if ( m[1] == *p ) {
649  *p = *m;
650  *m = *--stop;
651  m[1] = *--stop;
652  break;
653  }
654  else m += 2;
655  } while ( m < stop );
656  }
657  p++;
658  }
659  return(retval);
660  }
661  if ( number <= 2 ) return(0);
662  else { goto NextE4; }
663 }
664 
665 /*
666  #] Trace4no :
667  #[ Trace4 : WORD Trace4(term,params,num,level)
668 
669  Generates traces of the string of gamma matrices in 'instring'.
670 
671  The difference with the routine tracen ( for n dimensions )
672  lies in the treatment of gamma 5 and the specific form
673  of the Chisholm identities. The identities used here are
674  g(mu)*g(a1)*...*g(an)*g(mu)=
675  n=odd: -2*g(an)*...*g(a1) ( reversed order )
676  n=even: 2*g(an)*g(a1)*...*g(a(n-1))
677  +2*g(a(n-1))*...*g(a1)*g(an)
678  There is a special case for n=2 : 4*d(a1,a2)*gi
679 
680  The main difference with the old fortran version lies in
681  the recursion that is used here. That cleans up the variables
682  very much.
683 
684  The contents of the AT.TMout array are:
685  length,type,gamma5,gamma's
686 
687  The space for the vectors in t is at most 14 * number.
688 
689  The condition params[5] == 0 corresponds to finding gamma6*gamma7
690  during the pick up of the matrices.
691 
692 */
693 
694 WORD Trace4(PHEAD WORD *term, WORD *params, WORD num, WORD level)
695 {
696  GETBIDENTITY
697  TRACES *t;
698  WORD *p, *m, number, i;
699  WORD *OldW;
700  WORD j, minimum, minimum2, *min, *stopper;
701  OldW = AT.WorkPointer;
702  if ( AN.numtracesctack >= AN.intracestack ) {
703  number = AN.intracestack + 2;
704  t = (TRACES *)Malloc1(number*sizeof(TRACES),"TRACES-struct");
705  if ( AN.tracestack ) {
706  for ( i = 0; i < AN.intracestack; i++ ) { t[i] = AN.tracestack[i]; }
707  M_free(AN.tracestack,"TRACES-struct");
708  }
709  AN.tracestack = t;
710  AN.intracestack = number;
711  }
712 
713  number = *params - 6;
714  if ( number < 0 || ( number & 1 ) || !params[5] ) return(0);
715 
716  t = AN.tracestack + AN.numtracesctack;
717  AN.numtracesctack++;
718 
719  t->finalstep = ( params[2] & 16 ) ? 1 : 0;
720  t->gamma5 = params[3];
721  if ( t->finalstep && t->gamma5 != GAMMA1 ) {
722  MLOCK(ErrorMessageLock);
723  MesPrint("Gamma5 not allowed in this option of the trace command");
724  MUNLOCK(ErrorMessageLock);
725  AN.numtracesctack--;
726  SETERROR(-1)
727  }
728  t->inlist = AT.WorkPointer;
729  t->accup = t->accu = t->inlist + number;
730  t->perm = t->accu + (number*2);
731  t->eers = t->perm + number;
732  if ( ( AT.WorkPointer += 19 * number ) >= AT.WorkTop ) {
733  MLOCK(ErrorMessageLock);
734  MesWork();
735  MUNLOCK(ErrorMessageLock);
736  return(-1);
737  }
738  t->num = num;
739  t->level = level;
740  p = t->inlist;
741  m = params+6;
742  for ( i = 0; i < number; i++ ) *p++ = *m++;
743  t->termp = term;
744  t->factor = params[4];
745  t->allsign = params[5];
746  if ( number >= 10 || ( t->gamma5 != GAMMA1 && number > 4 ) ) {
747 /*
748  The next code should `normal order' the string
749  We need the lexicographic smallest string, taking also the
750  reverse strings into account.
751 */
752  minimum = 0; min = t->inlist;
753  stopper = min + number;
754  for ( i = 1; i < number; i++ ) {
755  p = min;
756  m = t->inlist + i;
757  for ( j = 0; j < number; j++ ) {
758  if ( *p < *m ) break;
759  if ( *p > *m ) {
760  min = t->inlist+i;
761  minimum = i;
762  break;
763  }
764  p++; m++;
765  if ( m >= stopper ) m = t->inlist;
766  if ( p >= stopper ) p = t->inlist;
767  }
768  }
769  p = min;
770  min = m = AT.WorkPointer;
771  i = number;
772  while ( --i >= 0 ) {
773  *m++ = *p++;
774  if ( p >= stopper ) p = t->inlist;
775  }
776  p = t->inlist;
777  m = min;
778  i = number;
779  while ( --i >= 0 ) *p++ = *m++;
780  p = t->inlist;
781  m = stopper;
782  while ( p < m ) { /* reverse string */
783  i = *p; *p++ = *--m; *m = i;
784  }
785  minimum2 = 0;
786  for ( i = 0; i < number; i++ ) {
787  p = min;
788  m = t->inlist + i;
789  for ( j = 0; j < number; j++ ) {
790  if ( *p < *m ) break;
791  if ( *p > *m ) {
792  m = t->inlist + i;
793  p = min;
794  j = number;
795  while ( --j >= 0 ) {
796  *p++ = *m++;
797  if ( m >= stopper ) m = t->inlist;
798  }
799  minimum2 = i;
800  break;
801  }
802  p++; m++;
803  if ( m >= stopper ) m = t->inlist;
804  }
805  }
806  minimum ^= minimum2;
807  if ( ( minimum & 1 ) != 0 ) {
808  if ( t->gamma5 == GAMMA5 ) t->allsign = - t->allsign;
809  else if ( t->gamma5 != GAMMA1 )
810  t->gamma5 = GAMMA6 + GAMMA7 - t->gamma5;
811  }
812  p = min; m = t->inlist; i = number;
813  while ( --i >= 0 ) *m++ = *p++;
814 /*
815  Now the trace is in normal order
816 */
817  }
818  number = Trace4Gen(BHEAD t,number);
819  AT.WorkPointer = OldW;
820  AN.numtracesctack--;
821  return(number);
822 }
823 
824 /*
825  #] Trace4 :
826  #[ Trace4Gen : WORD Trace4Gen(t,number)
827 
828  The recursive breakdown of a trace in 4 dimensions.
829  We test first whether the trace has zero or two gamma's left.
830  This case can be done quickly.
831  Next we test whether we can eliminate adjacent objects that are
832  the same.
833  Then we test for Chisholm identities (I). First for identities
834  with an odd number of gamma matrices (II), then for those with an
835  even number of matrices (III). The special thing here is the demand
836  that the contraction be between indices with 4 dimensions only.
837  Then there is a scan for objects that are the same, not regarding
838  their type (IV). This is exactly the same as in n dimensions.
839  Finally we have a string left in which all objects are different (V).
840  This case is treated by the routine Trace4no (no stands for no
841  objects are the same).
842 
843  In case I we have one d_ of which the result of the contraction
844  has not yet been fixed.
845  Case II gives just a reordering of the matrices and a factor -2.
846  Case III gives two terms: one for the anti commutation, such that
847  the number of intermediate matrices becomes odd and the other
848  from the Chisholm rule for an odd number of matrices. Both have
849  a factor 2.
850  Case IV gives m+1 terms when m is the number of matrices inbetween.
851  We take the shortest path. The sign alternates and all terms have
852  a factor two, except for the last one.
853 
854 */
855 
856 WORD Trace4Gen(PHEAD TRACES *t, WORD number)
857 {
858  GETBIDENTITY
859  WORD *termout, *stop;
860  WORD *p, *m, oldval;
861  WORD *pold, *mold, diff, *oldstring, cp;
862 /*
863  #[ Special cases :
864 */
865  if ( number <= 2 ) { /* Special cases */
866  if ( t->gamma5 == GAMMA5 ) return(0);
867  termout = AT.WorkPointer;
868  p = t->termp;
869  stop = p + *p;
870  m = termout;
871  p++;
872  if ( p < stop ) do {
873  if ( *p == SUBEXPRESSION && p[2] == t->num ) {
874  oldstring = p;
875  p = t->termp;
876  do { *m++ = *p++; } while ( p < oldstring );
877  p += p[1];
878  *m++ = AC.lUniTrace[0];
879  *m++ = AC.lUniTrace[1];
880  *m++ = AC.lUniTrace[2];
881  *m++ = AC.lUniTrace[3];
882  if ( number == 2 || t->accup > t->accu ) {
883  oldstring = m;
884  *m++ = DELTA;
885  *m++ = 4;
886  if ( number == 2 ) {
887  *m++ = t->inlist[0];
888  *m++ = t->inlist[1];
889  }
890  if ( t->accup > t->accu ) {
891  pold = p;
892  p = t->accu;
893  while ( p < t->accup ) *m++ = *p++;
894  oldstring[1] = WORDDIF(m,oldstring);
895  p = pold;
896  }
897  }
898  if ( t->factor ) {
899  *m++ = SNUMBER;
900  *m++ = 4;
901  *m++ = 2;
902  *m++ = t->factor;
903  }
904  do { *m++ = *p++; } while ( p < stop );
905  *termout = WORDDIF(m,termout);
906  if ( t->allsign < 0 ) m[-1] = -m[-1];
907  if ( ( AT.WorkPointer = m ) > AT.WorkTop ) {
908  MLOCK(ErrorMessageLock);
909  MesWork();
910  MUNLOCK(ErrorMessageLock);
911  return(-1);
912  }
913  *AN.RepPoint = 1;
914  AR.expchanged = 1;
915  if ( *termout ) {
916  *AN.RepPoint = 1;
917  AR.expchanged = 1;
918  if ( Generator(BHEAD termout,t->level) ) goto TracCall;
919  }
920  AT.WorkPointer= termout;
921  return(0);
922  }
923  p += p[1];
924  } while ( p < stop );
925  return(0);
926  }
927 /*
928  #] Special cases :
929  #[ Adjacent objects :
930 */
931  p = t->inlist;
932  stop = p + number - 1;
933  if ( *p == *stop ) { /* First and last of string */
934  oldval = *p;
935  *(t->accup)++ = *p;
936  *(t->accup)++ = *p;
937  m = p+1;
938  while ( m < stop ) *p++ = *m++;
939  if ( t->gamma5 != GAMMA1 ) {
940  if ( t->gamma5 == GAMMA5 ) t->allsign = - t->allsign;
941  else if ( t->gamma5 == GAMMA6 ) t->gamma5 = GAMMA7;
942  else if ( t->gamma5 == GAMMA7 ) t->gamma5 = GAMMA6;
943  }
944  if ( Trace4Gen(BHEAD t,number-2) ) goto TracCall;
945  t = AN.tracestack + AN.numtracesctack - 1;
946  if ( t->gamma5 != GAMMA1 ) {
947  if ( t->gamma5 == GAMMA5 ) t->allsign = - t->allsign;
948  else if ( t->gamma5 == GAMMA6 ) t->gamma5 = GAMMA7;
949  else if ( t->gamma5 == GAMMA7 ) t->gamma5 = GAMMA6;
950  }
951  while ( p > t->inlist ) *--m = *--p;
952  *p = *stop = oldval;
953  t->accup -= 2;
954  return(0);
955  }
956  do {
957  if ( *p == p[1] ) { /* Adjacent in string */
958  oldval = *p;
959  pold = p;
960  m = p+2;
961  *(t->accup)++ = *p;
962  *(t->accup)++ = *p;
963  while ( m <= stop ) *p++ = *m++;
964  if ( Trace4Gen(BHEAD t,number-2) ) goto TracCall;
965  t = AN.tracestack + AN.numtracesctack - 1;
966  while ( p > pold ) *--m = *--p;
967  *p++ = oldval;
968  *p++ = oldval;
969  t->accup -= 2;
970  return(0);
971  }
972  p++;
973  } while ( p < stop );
974 /*
975  #] Adjacent objects :
976  #[ Odd Contraction :
977 */
978  p = t->inlist;
979  stop = p + number;
980  do {
981  if ( *p >= AM.OffsetIndex && (
982  ( *p < WILDOFFSET + AM.OffsetIndex &&
983  indices[*p-AM.OffsetIndex].dimension == 4 )
984  || ( *p >= WILDOFFSET + AM.OffsetIndex && AC.lDefDim == 4 ) ) ) {
985  m = p+2;
986  while ( m < stop ) {
987  if ( *p == *m ) {
988  pold = p;
989  mold = m;
990  oldval = *p;
991  (t->factor)++;
992  t->allsign = - t->allsign;
993  *p++ = *--m;
994  m--;
995  while ( m > p ) { diff = *p; *p++ = *m; *m-- = diff; }
996  p = mold - 1;
997  m = mold + 1;
998  while ( m < stop ) *p++ = *m++;
999  if ( Trace4Gen(BHEAD t,number-2) ) goto TracCall;
1000  t = AN.tracestack + AN.numtracesctack - 1;
1001  m--;
1002  while ( m > mold ) *m-- = *--p;
1003  p = pold;
1004  *m-- = oldval;
1005  *m-- = *p;
1006  *p++ = oldval;
1007  while ( m > p ) { diff = *p; *p++ = *m; *m-- = diff; }
1008  t->allsign = - t->allsign;
1009  (t->factor)--;
1010  return(0);
1011  }
1012  m += 2;
1013  }
1014  }
1015  p++;
1016  } while ( p < stop );
1017 /*
1018  #] Odd Contraction :
1019  #[ Even Contraction :
1020  First the case with two matrices inbetween.
1021 */
1022  p = t->inlist;
1023  stop = p + number;
1024  do {
1025  if ( *p >= AM.OffsetIndex && (
1026  ( *p < WILDOFFSET + AM.OffsetIndex &&
1027  indices[*p-AM.OffsetIndex].dimension == 4 )
1028  || ( *p >= WILDOFFSET + AM.OffsetIndex && AC.lDefDim == 4 ) ) ) {
1029  m = p+3;
1030  if ( m >= stop ) m -= number;
1031  if ( *p == *m ) {
1032  WORD oldfactor, old5;
1033  oldstring = AT.WorkPointer;
1034  AT.WorkPointer += number;
1035  oldfactor = t->allsign;
1036  old5 = t->gamma5;
1037  if ( m < p ) cp = (WORDDIF(m,t->inlist) + 1 ) & 1;
1038  else cp = 0;
1039  if ( cp && ( t->gamma5 != GAMMA1 ) ) {
1040  if ( t->gamma5 == GAMMA5 ) t->allsign = -t->allsign;
1041  else if ( t->gamma5 == GAMMA6 ) t->gamma5 = GAMMA7;
1042  else if ( t->gamma5 == GAMMA7 ) t->gamma5 = GAMMA6;
1043  }
1044  mold = m;
1045  p = oldstring;
1046  m = t->inlist;
1047  while ( m < stop ) *p++ = *m++;
1048 /*
1049  Rotate the string
1050 */
1051  m = mold + 1;
1052  p = t->inlist;
1053  while ( m < stop ) *p++ = *m++;
1054  m = oldstring;
1055  if ( !cp && ((WORDDIF(stop,p))&1) != 0 && ( t->gamma5 != GAMMA1 ) ) {
1056  if ( t->gamma5 == GAMMA5 ) t->allsign = -t->allsign;
1057  else if ( t->gamma5 == GAMMA6 ) t->gamma5 = GAMMA7;
1058  else if ( t->gamma5 == GAMMA7 ) t->gamma5 = GAMMA6;
1059  }
1060  while ( p < stop ) *p++ = *m++;
1061  t->factor += 2;
1062  m = p - 3;
1063  p = t->inlist;
1064  oldval = number - 4;
1065  while ( oldval > 0 ) {
1066  if ( *p >= AM.OffsetIndex && (
1067  ( *p < WILDOFFSET + AM.OffsetIndex &&
1068  indices[*p-AM.OffsetIndex].dimension )
1069  || ( *p >= WILDOFFSET + AM.OffsetIndex && AC.lDefDim ) ) ) {
1070  if ( *p == *m ) {
1071  *p = m[1];
1072  break;
1073  }
1074  else if ( *p == m[1] ) {
1075  *p = *m;
1076  break;
1077  }
1078  }
1079  p++; oldval--;
1080  }
1081  if ( oldval <= 0 ) {
1082  *(t->accup)++ = *m++;
1083  *(t->accup)++ = *m++;
1084  }
1085  if ( Trace4Gen(BHEAD t,number-4) ) goto TracCall;
1086  t = AN.tracestack + AN.numtracesctack - 1;
1087  t->factor -= 2;
1088  if ( oldval <= 0 ) t->accup -= 2;
1089  t->gamma5 = old5;
1090  t->allsign = oldfactor;
1091  AT.WorkPointer = p = oldstring;
1092  m = t->inlist;
1093  while ( m < stop ) *m++ = *p++;
1094  return(0);
1095  }
1096  }
1097  p++;
1098  } while ( p < stop );
1099 /*
1100  The case with at least 4 matrices inbetween.
1101 */
1102  p = t->inlist;
1103  stop = p + number;
1104  do {
1105  if ( *p >= AM.OffsetIndex && (
1106  ( *p < WILDOFFSET + AM.OffsetIndex &&
1107  indices[*p-AM.OffsetIndex].dimension == 4 )
1108  || ( *p >= WILDOFFSET + AM.OffsetIndex && AC.lDefDim == 4 ) ) ) {
1109  m = p+5;
1110  while ( m < stop ) {
1111  if ( *p == *m ) {
1112  WORD *pex, *mex;
1113  pold = p;
1114  mold = m;
1115  oldval = *p;
1116 /*
1117  g_(1,mu)*g_(1,a1)*...*g_(1,aj)*g_(1,an)*g_(1,mu) ->
1118  first:
1119  2*g_(1,an)*g_(1,a1)*...*g_(1,aj)
1120 */
1121  (t->factor)++;
1122 /*
1123  The variable hulp seems unnecessary
1124  *p = hulp = m[-1];
1125 */
1126  *p = m[-1];
1127  p = m - 1;
1128  m++;
1129  while ( m < stop ) *p++ = *m++;
1130  if ( Trace4Gen(BHEAD t,number-2) ) goto TracCall;
1131  t = AN.tracestack + AN.numtracesctack - 1;
1132  pex = p; mex = m;
1133  p = pold;
1134  m = mold - 2;
1135  while ( m > p ) { diff = *p; *p++ = *m; *m-- = diff; }
1136 /*
1137  and then:
1138  2*g_(1,aj)*...*g_(1,a1)*g_(1,an)
1139 */
1140  if ( Trace4Gen(BHEAD t,number-2) ) goto TracCall;
1141  t = AN.tracestack + AN.numtracesctack - 1;
1142  p = pold;
1143  m = mold - 2;
1144  while ( m > p ) { diff = *p; *p++ = *m; *m-- = diff; }
1145  m = mex;
1146  p = pex;
1147  m--;
1148  while ( m > mold ) *m-- = *--p;
1149  m = mold;
1150  *m-- = oldval;
1151  p = pold;
1152  *m = *p;
1153  *p = oldval;
1154  (t->factor)--;
1155  return(0);
1156  }
1157  m += 2;
1158  }
1159  }
1160  p++;
1161  } while ( p < stop );
1162 /*
1163  #] Even Contraction :
1164  #[ Same Objects :
1165 */
1166  p = t->inlist;
1167  stop = p + number - 1;
1168  diff = 2;
1169  do {
1170  p = t->inlist;
1171  while ( p <= stop ) {
1172  m = p + diff;
1173  if ( m > stop ) m -= number;
1174  if ( *p == *m ) {
1175  WORD oldfactor, c, old5;
1176  oldfactor = t->allsign;
1177  old5 = t->gamma5;
1178  cp = (WORDDIF(m,t->inlist)) & 1;
1179  if ( !cp && ( t->gamma5 != GAMMA1 ) ) {
1180  if ( t->gamma5 == GAMMA5 ) t->allsign = -t->allsign;
1181  else if ( t->gamma5 == GAMMA6 ) t->gamma5 = GAMMA7;
1182  else if ( t->gamma5 == GAMMA7 ) t->gamma5 = GAMMA6;
1183  }
1184  oldstring = AT.WorkPointer;
1185  AT.WorkPointer += number;
1186  mold = m;
1187  oldval = *p;
1188  p = oldstring;
1189  m = t->inlist;
1190  while ( m <= stop ) *p++ = *m++;
1191 /*
1192  Rotate the string
1193 */
1194  m = mold + 1;
1195  p = t->inlist;
1196  while ( m <= stop ) *p++ = *m++;
1197  m = oldstring;
1198  while ( p <= stop ) *p++ = *m++;
1199  (t->factor)++;
1200  p -= diff + 1;
1201  m = stop;
1202  *(t->accup) = oldval;
1203  t->accup += 2;
1204  m--;
1205  while ( m > p ) {
1206  c = t->accup[-1];
1207  t->accup[-1] = *m;
1208  *m = c;
1209  if ( Trace4Gen(BHEAD t,number-2) ) goto Trac4Call;
1210  t = AN.tracestack + AN.numtracesctack - 1;
1211  m--;
1212  t->allsign = - t->allsign;
1213  }
1214  c = t->accup[-1];
1215  t->accup[-1] = *m;
1216  *m = c;
1217  (t->factor)--;
1218  if ( Trace4Gen(BHEAD t,number-2) ) goto Trac4Call;
1219  t = AN.tracestack + AN.numtracesctack - 1;
1220  t->accup -= 2;
1221  t->allsign = oldfactor;
1222  AT.WorkPointer = p = oldstring;
1223  m = t->inlist;
1224  while ( m <= stop ) *m++ = *p++;
1225  t->gamma5 = old5;
1226  return(0);
1227  }
1228  p++;
1229  }
1230  } while ( ++diff <= (number>>1) );
1231 /*
1232  #] Same Objects :
1233  #[ All Different :
1234 
1235  Here we have a string with all different objects.
1236 
1237 */
1238  t->sgn = 0;
1239  termout = AT.WorkPointer;
1240  for(;;) {
1241  if ( t->finalstep == 0 ) diff = Trace4no(number,t->accup,t);
1242  else diff = TraceNno(number,t->accup,t);
1243 /* while ( ( diff = Trace4no(number,t->accup,t) ) != 0 ) */
1244  if ( diff == 0 ) break;
1245  p = t->termp;
1246  stop = p + *p;
1247  m = termout;
1248  p++;
1249  if ( p < stop ) do {
1250  if ( *p == SUBEXPRESSION && p[2] == t->num ) {
1251  oldstring = p;
1252  p = t->termp;
1253  do { *m++ = *p++; } while ( p < oldstring );
1254  p += p[1];
1255  pold = p;
1256  *m++ = AC.lUniTrace[0];
1257  *m++ = AC.lUniTrace[1];
1258  *m++ = AC.lUniTrace[2];
1259  *m++ = AC.lUniTrace[3];
1260  *m++ = SNUMBER;
1261  *m++ = 4;
1262  *m++ = 2;
1263  *m++ = t->factor;
1264  p = t->accup;
1265  oldval = number;
1266  if ( diff == 2 || diff == -2 ) {
1267  *m++ = LEVICIVITA;
1268  *m++ = 4+FUNHEAD;
1269  *m++ = DIRTYFLAG;
1270  FILLFUN3(m)
1271  *m++ = *p++; *m++ = *p++; *m++ = *p++; *m++ = *p++;
1272  oldval -= 4;
1273  }
1274  if ( oldval > 0 || t->accup > t->accu ) {
1275  oldstring = m;
1276  *m++ = DELTA;
1277  *m++ = oldval + 2;
1278  if ( oldval > 0 ) NCOPY(m,p,oldval);
1279  if ( t->accup > t->accu ) {
1280  p = t->accu;
1281  while ( p < t->accup ) *m++ = *p++;
1282  oldstring[1] = WORDDIF(m,oldstring);
1283  }
1284  }
1285  p = pold;
1286  do { *m++ = *p++; } while ( p < stop );
1287  *termout = WORDDIF(m,termout);
1288  if ( ( diff ^ t->allsign ) < 0 ) m[-1] = - m[-1];
1289  if ( ( AT.WorkPointer = m ) > AT.WorkTop ) {
1290  MLOCK(ErrorMessageLock);
1291  MesWork();
1292  MUNLOCK(ErrorMessageLock);
1293  return(-1);
1294  }
1295  if ( *termout ) {
1296  *AN.RepPoint = 1;
1297  AR.expchanged = 1;
1298  if ( Generator(BHEAD termout,t->level) ) {
1299  AT.WorkPointer = termout;
1300  goto TracCall;
1301  }
1302  t = AN.tracestack + AN.numtracesctack - 1;
1303  }
1304  break;
1305  }
1306  p += p[1];
1307  } while ( p < stop );
1308  }
1309  AT.WorkPointer = termout;
1310  return(0);
1311 
1312 /*
1313  #] All Different :
1314 */
1315 Trac4Call:
1316  AT.WorkPointer = oldstring;
1317 TracCall:
1318  if ( AM.tracebackflag ) {
1319  MLOCK(ErrorMessageLock);
1320  MesCall("Trace4Gen");
1321  MUNLOCK(ErrorMessageLock);
1322  }
1323  return(-1);
1324 }
1325 
1326 /*
1327  #] Trace4Gen :
1328  #[ TraceNno : WORD TraceNno(number,kron,t)
1329 
1330  Routine takes the trace in N dimensions of a string
1331  of gamma matrices. It is assumed that there are no
1332  contractions and no vectors that are the same. For
1333  the treatment of those cases there are special routines,
1334  that call this routine as a final stage.
1335  The calling routine must reserve 'number' WORDs for perm
1336  and kron each.
1337  kron and perm may not be altered during the generation!
1338 
1339 */
1340 
1341 WORD TraceNno(WORD number, WORD *kron, TRACES *t)
1342 {
1343  WORD i, j, a, *p;
1344  if ( !t->sgn ) {
1345  if ( !number || ( number & 1 ) ) return(0);
1346  p = t->inlist;
1347  for ( i = 0; i < number; i++ ) {
1348  t->perm[i] = i;
1349  *kron++ = *p++;
1350  }
1351  t->sgn = 1;
1352  return(1);
1353  }
1354  else {
1355  i = number - 3;
1356  while ( i > 0 ) {
1357  a = kron[i];
1358  p = t->perm + i;
1359  for ( j = i + 1; j <= *p; j++ ) kron[j-1] = kron[j];
1360  kron[(*p)++] = a;
1361  if ( *p < number ) {
1362  a = kron[*p];
1363  j = *p;
1364  while ( j >= (i+1) ) { kron[j] = kron[j-1]; j--; }
1365  kron[i] = a;
1366  number -= 2;
1367  for ( j = i+2; j < number; j += 2 ) t->perm[j] = j;
1368  t->sgn = - t->sgn;
1369  return(t->sgn);
1370  }
1371  i -= 2;
1372  }
1373  }
1374  return(0);
1375 }
1376 
1377 /*
1378  #] TraceNno :
1379  #[ TraceN : WORD TraceN(term,params,num,level)
1380 */
1381 
1382 WORD TraceN(PHEAD WORD *term, WORD *params, WORD num, WORD level)
1383 {
1384  GETBIDENTITY
1385  TRACES *t;
1386  WORD *p, *m, number, i;
1387  WORD *OldW;
1388  if ( params[3] != GAMMA1 ) {
1389  MLOCK(ErrorMessageLock);
1390  MesPrint("Gamma5 not allowed in n-trace");
1391  MUNLOCK(ErrorMessageLock);
1392  SETERROR(-1)
1393  }
1394  OldW = AT.WorkPointer;
1395  if ( AN.numtracesctack >= AN.intracestack ) {
1396  number = AN.intracestack + 2;
1397  t = (TRACES *)Malloc1(number*sizeof(TRACES),"TRACES-struct");
1398  if ( AN.tracestack ) {
1399  for ( i = 0; i < AN.intracestack; i++ ) { t[i] = AN.tracestack[i]; }
1400  M_free(AN.tracestack,"TRACES-struct");
1401  }
1402  AN.tracestack = t;
1403  AN.intracestack = number;
1404  }
1405  number = *params - 6;
1406  if ( number < 0 || ( number & 1 ) || !params[5] ) return(0);
1407 
1408  t = AN.tracestack + AN.numtracesctack;
1409  AN.numtracesctack++;
1410 
1411  t->inlist = AT.WorkPointer;
1412  t->accup = t->accu = t->inlist + number;
1413  t->perm = t->accu + number;
1414  if ( ( AT.WorkPointer += 3 * number ) >= AT.WorkTop ) {
1415  AN.numtracesctack--;
1416  MLOCK(ErrorMessageLock);
1417  MesWork();
1418  MUNLOCK(ErrorMessageLock);
1419  return(-1);
1420  }
1421  t->num = num;
1422  t->level = level;
1423  p = t->inlist;
1424  m = params+6;
1425  for ( i = 0; i < number; i++ ) *p++ = *m++;
1426  t->termp = term;
1427  t->factor = params[4];
1428  t->allsign = params[5];
1429  number = TraceNgen(BHEAD t,number);
1430  AT.WorkPointer = OldW;
1431  AN.numtracesctack--;
1432  return(number);
1433 }
1434 
1435 /*
1436  #] TraceN :
1437  #[ TraceNgen : WORD TraceNgen(t,number)
1438 
1439  This routine is a simplified version of Trace4Gen. We know here
1440  only three cases: Adjacent objects, same objects and all different.
1441  The othere difference lies of course in the struct which is now
1442  not of type TRACES, but of type TRACES.
1443 
1444 */
1445 
1446 WORD TraceNgen(PHEAD TRACES *t, WORD number)
1447 {
1448  GETBIDENTITY
1449  WORD *termout, *stop;
1450  WORD *p, *m, oldval;
1451  WORD *pold, *mold, diff, *oldstring;
1452 /*
1453  #[ Special cases :
1454 */
1455  if ( number <= 2 ) { /* Special cases */
1456  termout = AT.WorkPointer;
1457  p = t->termp;
1458  stop = p + *p;
1459  m = termout;
1460  p++;
1461  if ( p < stop ) do {
1462  if ( *p == SUBEXPRESSION && p[2] == t->num ) {
1463  oldstring = p;
1464  p = t->termp;
1465  do { *m++ = *p++; } while ( p < oldstring );
1466  p += p[1];
1467  *m++ = AC.lUniTrace[0];
1468  *m++ = AC.lUniTrace[1];
1469  *m++ = AC.lUniTrace[2];
1470  *m++ = AC.lUniTrace[3];
1471  if ( number == 2 || t->accup > t->accu ) {
1472  oldstring = m;
1473  *m++ = DELTA;
1474  *m++ = 4;
1475  if ( number == 2 ) {
1476  *m++ = t->inlist[0];
1477  *m++ = t->inlist[1];
1478  }
1479  if ( t->accup > t->accu ) {
1480  pold = p;
1481  p = t->accu;
1482  while ( p < t->accup ) *m++ = *p++;
1483  oldstring[1] = WORDDIF(m,oldstring);
1484  p = pold;
1485  }
1486  }
1487  if ( t->factor ) {
1488  *m++ = SNUMBER;
1489  *m++ = 4;
1490  *m++ = 2;
1491  *m++ = t->factor;
1492  }
1493  do { *m++ = *p++; } while ( p < stop );
1494  *termout = WORDDIF(m,termout);
1495  if ( t->allsign < 0 ) m[-1] = -m[-1];
1496  if ( ( AT.WorkPointer = m ) > AT.WorkTop ) {
1497  MLOCK(ErrorMessageLock);
1498  MesWork();
1499  MUNLOCK(ErrorMessageLock);
1500  return(-1);
1501  }
1502  if ( *termout ) {
1503  *AN.RepPoint = 1;
1504  AR.expchanged = 1;
1505  if ( Generator(BHEAD termout,t->level) ) goto TracCall;
1506  }
1507  AT.WorkPointer= termout;
1508  return(0);
1509  }
1510  p += p[1];
1511  } while ( p < stop );
1512  return(0);
1513  }
1514 /*
1515  #] Special cases :
1516  #[ Adjacent objects :
1517 */
1518  p = t->inlist;
1519  stop = p + number - 1;
1520  if ( *p == *stop ) { /* First and last of string */
1521  oldval = *p;
1522  *(t->accup)++ = *p;
1523  *(t->accup)++ = *p;
1524  m = p+1;
1525  while ( m < stop ) *p++ = *m++;
1526  if ( TraceNgen(BHEAD t,number-2) ) goto TracCall;
1527  t = AN.tracestack + AN.numtracesctack - 1;
1528  while ( p > t->inlist ) *--m = *--p;
1529  *p = *stop = oldval;
1530  t->accup -= 2;
1531  return(0);
1532  }
1533  do {
1534  if ( *p == p[1] ) { /* Adjacent in string */
1535  oldval = *p;
1536  pold = p;
1537  m = p+2;
1538  *(t->accup)++ = *p;
1539  *(t->accup)++ = *p;
1540  while ( m <= stop ) *p++ = *m++;
1541  if ( TraceNgen(BHEAD t,number-2) ) goto TracCall;
1542  t = AN.tracestack + AN.numtracesctack - 1;
1543  while ( p > pold ) *--m = *--p;
1544  *p++ = oldval;
1545  *p++ = oldval;
1546  t->accup -= 2;
1547  return(0);
1548  }
1549  p++;
1550  } while ( p < stop );
1551 /*
1552  #] Adjacent objects :
1553  #[ Same Objects :
1554 */
1555  p = t->inlist;
1556  stop = p + number - 1;
1557  diff = 2;
1558  do {
1559  p = t->inlist;
1560  while ( p <= stop ) {
1561  m = p + diff;
1562  if ( m > stop ) m -= number;
1563  if ( *p == *m ) {
1564  WORD oldfactor, c;
1565  oldstring = AT.WorkPointer;
1566  AT.WorkPointer += number;
1567  mold = m;
1568  oldval = *p;
1569  p = oldstring;
1570  m = t->inlist;
1571  while ( m <= stop ) *p++ = *m++;
1572 /*
1573  Rotate the string
1574 */
1575  {
1576  m = mold + 1;
1577  p = t->inlist;
1578  while ( m <= stop ) *p++ = *m++;
1579  m = oldstring;
1580  while ( p <= stop ) *p++ = *m++;
1581  }
1582  oldfactor = t->allsign;
1583  (t->factor)++;
1584  p -= diff + 1;
1585  m = stop;
1586  if ( oldval >= ( AM.OffsetIndex + WILDOFFSET ) ||
1587  ( oldval >= AM.OffsetIndex
1588  && indices[oldval-AM.OffsetIndex].dimension ) ) {
1589  m--;
1590 /*
1591  We distinguish 4 cases:
1592  m-p=1 Use g_(1,mu,a,mu) = (2-d_(mu,mu))*g_(1,a)
1593  m-p=2 Use g_(1,mu,a,b,mu) = 4*d_(a,b)+(d_(mu,mu)-4)*g_(1,a,b)
1594  m-p=3 Use g_(1,mu,a,b,c,mu) = -2*g_(1,c,b,a)
1595  -(d_(mu,mu)-4)*g_(1,a,b,c)
1596  m-p>3 Reduce down to m-p=3 with the old technique
1597 */
1598  while ( m > (p+3) ) {
1599  c = *p;
1600  *p = *m;
1601  *m = c;
1602  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1603  t = AN.tracestack + AN.numtracesctack - 1;
1604  m--;
1605  t->allsign = - t->allsign;
1606  }
1607  switch ( WORDDIF(m,p) ) {
1608  case 1:
1609  c = *p;
1610  *p = *m;
1611  *m = c;
1612  if ( oldval < ( AM.OffsetIndex + WILDOFFSET )
1613  && indices[oldval-AM.OffsetIndex].nmin4
1614  != -NMIN4SHIFT ) {
1615  t->allsign = - t->allsign;
1616  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1617  t = AN.tracestack + AN.numtracesctack - 1;
1618  (t->factor)--;
1619  *(t->accup)++ = SUMMEDIND;
1620  *(t->accup)++ =
1621  indices[oldval-AM.OffsetIndex].nmin4;
1622  }
1623  else
1624  {
1625  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1626  t = AN.tracestack + AN.numtracesctack - 1;
1627  t->allsign = - t->allsign;
1628  (t->factor)--;
1629  *(t->accup)++ = oldval;
1630  *(t->accup)++ = oldval;
1631  }
1632  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1633  t = AN.tracestack + AN.numtracesctack - 1;
1634  t->accup -= 2;
1635  break;
1636  case 2:
1637  { WORD one, two;
1638  one = *p = p[1];
1639  two = p[1] = *m;
1640  (t->factor)++; /* 4 */
1641  *(t->accup)++ = *p; /* d_(a,b) */
1642  *(t->accup)++ = *m;
1643  if ( TraceNgen(BHEAD t,number-4) ) goto TracnCall;
1644  t = AN.tracestack + AN.numtracesctack - 1;
1645  *p = one; p[1] = two;
1646  t->accup -= 2;
1647  if ( oldval < ( AM.OffsetIndex + WILDOFFSET )
1648  && indices[oldval-AM.OffsetIndex].nmin4
1649  != -NMIN4SHIFT ) {
1650  t->factor -= 2;
1651  *(t->accup)++ = SUMMEDIND;
1652  *(t->accup)++ =
1653  indices[oldval-AM.OffsetIndex].nmin4;
1654  }
1655  else {
1656  t->allsign = - t->allsign;
1657  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1658  t = AN.tracestack + AN.numtracesctack - 1;
1659  t->allsign = - t->allsign;
1660  t->factor -= 2;
1661  *(t->accup)++ = oldval;
1662  *(t->accup)++ = oldval;
1663  }
1664  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1665  t = AN.tracestack + AN.numtracesctack - 1;
1666  t->accup -= 2;
1667  }
1668  break;
1669  default:
1670  c = *p;
1671  *p = *m;
1672  *m = c;
1673  c = m[-1]; m[-1] = m[-2]; m[-2] = c;
1674  t->allsign = - t->allsign;
1675  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1676  t = AN.tracestack + AN.numtracesctack - 1;
1677  m--;
1678  c = *p;
1679  *p = *m;
1680  *m = c;
1681  (t->factor)--;
1682  if ( oldval < ( AM.OffsetIndex + WILDOFFSET )
1683  && indices[oldval-AM.OffsetIndex].nmin4
1684  != -NMIN4SHIFT ) {
1685  *(t->accup)++ = SUMMEDIND;
1686  *(t->accup)++ =
1687  indices[oldval-AM.OffsetIndex].nmin4;
1688  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1689  t = AN.tracestack + AN.numtracesctack - 1;
1690  t->accup -= 2;
1691  t->allsign = - t->allsign;
1692  }
1693  else
1694  {
1695  *(t->accup)++ = oldval;
1696  *(t->accup)++ = oldval;
1697  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1698  t = AN.tracestack + AN.numtracesctack - 1;
1699  t->accup -= 2;
1700  t->allsign = - t->allsign;
1701  t->factor += 2;
1702  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1703  t = AN.tracestack + AN.numtracesctack - 1;
1704  t->factor -= 2;
1705  }
1706  break;
1707  }
1708  }
1709  else {
1710  *(t->accup) = oldval;
1711  t->accup += 2;
1712  m--;
1713  while ( m > p ) {
1714  c = t->accup[-1];
1715  t->accup[-1] = *m;
1716  *m = c;
1717  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1718  t = AN.tracestack + AN.numtracesctack - 1;
1719  m--;
1720  t->allsign = - t->allsign;
1721  }
1722  c = t->accup[-1];
1723  t->accup[-1] = *m;
1724  *m = c;
1725  (t->factor)--;
1726  if ( TraceNgen(BHEAD t,number-2) ) goto TracnCall;
1727  t = AN.tracestack + AN.numtracesctack - 1;
1728  t->accup -= 2;
1729  }
1730  t->allsign = oldfactor;
1731  p = oldstring;
1732  m = t->inlist;
1733  while ( m <= stop ) *m++ = *p++;
1734  AT.WorkPointer = oldstring;
1735  return(0);
1736  }
1737  p++;
1738  }
1739  diff++;
1740  } while ( diff <= (number>>1) );
1741 /*
1742  #] Same Objects :
1743  #[ All Different :
1744 
1745  Here we have a string with all different objects.
1746 
1747 */
1748  t->sgn = 0;
1749  termout = AT.WorkPointer;
1750  while ( ( diff = TraceNno(number,t->accup,t) ) != 0 ) {
1751  p = t->termp;
1752  stop = p + *p;
1753  m = termout;
1754  p++;
1755  if ( p < stop ) do {
1756  if ( *p == SUBEXPRESSION && p[2] == t->num ) {
1757  oldstring = p;
1758  p = t->termp;
1759  do { *m++ = *p++; } while ( p < oldstring );
1760  p += p[1];
1761  pold = p;
1762  *m++ = AC.lUniTrace[0];
1763  *m++ = AC.lUniTrace[1];
1764  *m++ = AC.lUniTrace[2];
1765  *m++ = AC.lUniTrace[3];
1766  *m++ = SNUMBER;
1767  *m++ = 4;
1768  *m++ = 2;
1769  *m++ = t->factor;
1770  p = t->accup;
1771  oldval = number;
1772  oldstring = m;
1773  *m++ = DELTA;
1774  *m++ = oldval + 2;
1775  NCOPY(m,p,oldval);
1776  if ( t->accup > t->accu ) {
1777  p = t->accu;
1778  while ( p < t->accup ) *m++ = *p++;
1779  oldstring[1] = WORDDIF(m,oldstring);
1780  }
1781  p = pold;
1782  do { *m++ = *p++; } while ( p < stop );
1783  *termout = WORDDIF(m,termout);
1784  if ( ( diff ^ t->allsign ) < 0 ) m[-1] = - m[-1];
1785  if ( ( AT.WorkPointer = m ) > AT.WorkTop ) {
1786  MLOCK(ErrorMessageLock);
1787  MesWork();
1788  MUNLOCK(ErrorMessageLock);
1789  return(-1);
1790  }
1791  if ( *termout ) {
1792  *AN.RepPoint = 1;
1793  AR.expchanged = 1;
1794  if ( Generator(BHEAD termout,t->level) ) {
1795  AT.WorkPointer = termout;
1796  goto TracCall;
1797  }
1798  t = AN.tracestack + AN.numtracesctack - 1;
1799  }
1800  break;
1801  }
1802  p += p[1];
1803  } while ( p < stop );
1804  }
1805  AT.WorkPointer = termout;
1806  return(0);
1807 
1808 /*
1809  #] All Different :
1810 */
1811 TracnCall:
1812  AT.WorkPointer = oldstring;
1813 TracCall:
1814  if ( AM.tracebackflag ) {
1815  MLOCK(ErrorMessageLock);
1816  MesCall("TraceNGen");
1817  MUNLOCK(ErrorMessageLock);
1818  }
1819  return(-1);
1820 }
1821 
1822 /*
1823  #] TraceNgen :
1824  #[ Traces : WORD Traces(term,params,num,level)
1825 
1826  The contents of the AT.TMout array are:
1827  length,type,subtype,gamma5,factor,sign,gamma's
1828 
1829 */
1830 
1831 WORD Traces(PHEAD WORD *term, WORD *params, WORD num, WORD level)
1832 {
1833  GETBIDENTITY
1834  switch ( AT.TMout[2] ) { /* Subtype gives dimension */
1835  case 0:
1836  return(TraceN(BHEAD term,params,num,level));
1837  case 4:
1838  return(Trace4(BHEAD term,params,num,level));
1839  case 12:
1840  return(Trace4(BHEAD term,params,num,level));
1841  case 20:
1842  return(Trace4(BHEAD term,params,num,level));
1843  default:
1844  return(0);
1845  }
1846 }
1847 
1848 /*
1849  #] Traces :
1850  #[ TraceFind : WORD TraceFind(term,params)
1851 */
1852 
1853 WORD TraceFind(PHEAD WORD *term, WORD *params)
1854 {
1855  GETBIDENTITY
1856  WORD *p, *m, *to;
1857  WORD *termout, *stop, *stop2, number = 0;
1858  WORD first = 1;
1859  WORD type, spinline, sp;
1860  type = params[3];
1861  spinline = params[4];
1862  if ( spinline < 0 ) { /* $ variable. Evaluate */
1863  sp = DolToIndex(BHEAD -spinline);
1864  if ( AN.ErrorInDollar || sp < 0 ) {
1865  MLOCK(ErrorMessageLock);
1866  MesPrint("$%s does not have an index value in trace statement in module %l",
1867  DOLLARNAME(Dollars,-spinline),AC.CModule);
1868  MUNLOCK(ErrorMessageLock);
1869  return(0);
1870  }
1871  spinline = sp;
1872  }
1873  to = AT.TMout;
1874  to++;
1875  *to++ = TAKETRACE;
1876  *to++ = type;
1877  *to++ = GAMMA1;
1878  *to++ = 0; /* Powers of two */
1879  *to++ = 1; /* sign */
1880  p = term;
1881  m = p + *p - 1;
1882  stop = m - ABS(*m);
1883  termout = m = AT.WorkPointer;
1884  m++;
1885  p++;
1886  while ( p < stop ) {
1887  stop2 = p + p[1];
1888  if ( *p == GAMMA && p[FUNHEAD] == spinline ) {
1889  if ( first ) {
1890  *m++ = SUBEXPRESSION;
1891  *m++ = SUBEXPSIZE;
1892  *m++ = -1;
1893  *m++ = 1;
1894  *m++ = DUMMYBUFFER;
1895  FILLSUB(m)
1896  first = 0;
1897  }
1898  p += FUNHEAD+1;
1899  while ( p < stop2 ) {
1900  if ( *p == GAMMA5 ) {
1901  if ( AT.TMout[3] == GAMMA5 ) AT.TMout[3] = GAMMA1;
1902  else if ( AT.TMout[3] == GAMMA1 ) AT.TMout[3] = GAMMA5;
1903  else if ( AT.TMout[3] == GAMMA7 ) AT.TMout[5] = -AT.TMout[5];
1904  if ( number & 1 ) AT.TMout[5] = - AT.TMout[5];
1905  p++;
1906  }
1907  else if ( *p == GAMMA6 ) {
1908  if ( number & 1 ) goto F7;
1909 F6: if ( AT.TMout[3] == GAMMA6 ) (AT.TMout[4])++;
1910  else if ( AT.TMout[3] == GAMMA1 ) AT.TMout[3] = GAMMA6;
1911  else if ( AT.TMout[3] == GAMMA5 ) AT.TMout[3] = GAMMA6;
1912  else if ( AT.TMout[3] == GAMMA7 ) AT.TMout[5] = 0;
1913  p++;
1914  }
1915  else if ( *p == GAMMA7 ) {
1916  if ( number & 1 ) goto F6;
1917 F7: if ( AT.TMout[3] == GAMMA7 ) (AT.TMout[4])++;
1918  else if ( AT.TMout[3] == GAMMA1 ) AT.TMout[3] = GAMMA7;
1919  else if ( AT.TMout[3] == GAMMA5 ) {
1920  AT.TMout[3] = GAMMA7;
1921  AT.TMout[5] = -AT.TMout[5];
1922  }
1923  else if ( AT.TMout[3] == GAMMA6 ) AT.TMout[5] = 0;
1924  p++;
1925  }
1926  else {
1927  *to++ = *p++;
1928  number++;
1929  }
1930  }
1931  }
1932  else {
1933  while ( p < stop2 ) *m++ = *p++;
1934  }
1935  }
1936  if ( first ) return(0);
1937  AT.TMout[0] = WORDDIF(to,AT.TMout);
1938  to = term;
1939  to += *to;
1940  while ( p < to ) *m++ = *p++;
1941  *termout = WORDDIF(m,termout);
1942  to = term;
1943  p = termout;
1944  do { *to++ = *p++; } while ( p < m );
1945  AT.WorkPointer = term + *term;
1946  return(1);
1947 }
1948 
1949 /*
1950  #] TraceFind :
1951  #[ Chisholm : WORD Chisholm(term,level,num)
1952 
1953  Routines for reorganizing traces.
1954  The command
1955  Chisholm,1;
1956  will collect the gamma matrices in spinline 1 and see whether
1957  they have an index in common with another gamma matrix. If this
1958  is the case the identity
1959  g_(2,mu)*Tr[g_(1,mu)*S(2)] = S(2)+SR(2)
1960  is applied (SR is the reversed string).
1961 */
1962 
1963 WORD Chisholm(PHEAD WORD *term, WORD level)
1964 {
1965  GETBIDENTITY
1966  WORD *t, *r, *m, *s, *tt, *rr;
1967  WORD *mat, *matpoint, *termout, *rdo;
1968  CBUF *C = cbuf+AM.rbufnum;
1969  WORD i, j, num = C->lhs[level][2], gam5;
1970  WORD norm = 0, k, *matp;
1971 /*
1972  #[ Find : Find possible matrices
1973 */
1974  mat = matpoint = AT.WorkPointer;
1975  t = term;
1976  r = t + *t - 1; r -= ABS(*r);
1977  t++;
1978  i = 0;
1979  gam5 = GAMMA1;
1980  while ( t < r ) {
1981  if ( *t == GAMMA && t[FUNHEAD] == num ) {
1982  m = t + t[1];
1983  t += FUNHEAD+1;
1984  while ( t < m ) {
1985  if ( *t >= 0 || *t < MINSPEC ) i++;
1986  else {
1987  if ( gam5 == GAMMA1 ) gam5 = *t;
1988  else if ( gam5 == GAMMA5 ) {
1989  if ( *t == GAMMA5 ) gam5 = GAMMA1;
1990  else if ( *t != GAMMA1 ) gam5 = *t;
1991  }
1992  }
1993  *matpoint++ = *t++;
1994  }
1995  }
1996  else t += t[1];
1997  }
1998  if ( ( i & 1 ) != 0 ) return(0); /* odd trace */
1999 /*
2000  #] Find :
2001  #[ Test : Test for contracted index
2002 
2003  This code should be modified.
2004 
2005  We have to check for all possible matches if C->lhs[level][3] == 1
2006  and the trace contains a gamma5, gamma6 or gamma7.
2007  Then we normalize by the number of possible contractions (norm) and
2008  do all of them. This way the Levi-Civita tensors have a maximum
2009  chance of cancelling each other. This option is activated with
2010  `contract' and `symmetrize'. Defaults are that they are on, but
2011  they can be switched off with nocontract and nosymmetrize.
2012 */
2013  s = mat;
2014  while ( s < matpoint ) {
2015 /*
2016  if ( *s < AM.OffsetIndex || ( *s < ( AM.OffsetIndex + WILDOFFSET ) &&
2017  indices[*s-AM.OffsetIndex].dimension == 0 ) ) {
2018 */
2019  if ( *s < AM.OffsetIndex || ( *s < ( AM.OffsetIndex + WILDOFFSET ) &&
2020  indices[*s-AM.OffsetIndex].dimension != 4 )
2021  || ( ( AC.lDefDim != 4 ) && ( *s >= ( AM.OffsetIndex + WILDOFFSET ) ) ) ) {
2022  s++; continue;
2023  }
2024  t = term+1;
2025  while ( t < r ) {
2026  if ( *t == GAMMA && t[FUNHEAD] != num ) {
2027  m = t + t[1];
2028  t += FUNHEAD+1;
2029  while ( t < m ) {
2030  if ( *t == *s ) {
2031  norm++;
2032  }
2033  t++;
2034  }
2035  }
2036  else t += t[1];
2037  }
2038  s++;
2039  }
2040  if ( norm == 0 ) return(Generator(BHEAD term,level)); /* No Action */
2041 /*
2042  #] Test :
2043  #[ Do : Process the string
2044 
2045  tt: The subterm
2046  t: The matrix
2047  s: The matrix in the relevant string
2048 
2049  Cycle the string in mat so that s is at the end.
2050  Copy the part till the critical GAMMA.
2051  Copy inside the critical string, copy S, copy tail inside string.
2052  Important to remember where S is so that we can reverse it later.
2053  Add term UnitTrace/2/norm.
2054  Copy rest of term.
2055  Continue execution with S.
2056  Reverse S.
2057  Continue execution with SR.
2058 */
2059 
2060  if ( C->lhs[level][3] == 0 /* || gam5 == GAMMA1 */ ) norm = 1;
2061 
2062  matp = matpoint;
2063  for ( k = 0; k < norm; k++ ) {
2064  matpoint = matp;
2065  s = mat;
2066  while ( s < matpoint ) {
2067 /*
2068  if ( *s < AM.OffsetIndex || ( *s < ( AM.OffsetIndex + WILDOFFSET ) &&
2069  indices[*s-AM.OffsetIndex].dimension == 0 ) ) {
2070 */
2071  if ( *s < AM.OffsetIndex || ( *s < ( AM.OffsetIndex + WILDOFFSET ) &&
2072  indices[*s-AM.OffsetIndex].dimension != 4 ) ) {
2073  s++; continue;
2074  }
2075  t = term+1;
2076  while ( t < r ) {
2077  if ( *t == GAMMA && t[FUNHEAD] != num ) {
2078  tt = t;
2079  m = t + t[1];
2080  t += FUNHEAD+1;
2081  while ( t < m ) {
2082  if ( *t == *s ) {
2083  i = WORDDIF(t,tt);
2084  m = mat;
2085  while ( m <= s ) *matpoint++ = *m++;
2086  t = mat;
2087  while ( m < matpoint ) *t++ = *m++;
2088  termout = t;
2089  m = termout + 1;
2090  t = term + 1;
2091  while ( t < tt ) {
2092  if ( *t != GAMMA || t[FUNHEAD] != num ) {
2093  j = t[1];
2094  NCOPY(m,t,j);
2095  }
2096  else t += t[1];
2097  }
2098 
2099  tt += tt[1];
2100  rdo = m;
2101  j = i;
2102  while ( --j >= 0 ) *m++ = *t++;
2103  matpoint = m;
2104  s = mat;
2105  while ( s < termout ) *m++ = *s++;
2106  m--;
2107  t++;
2108  while ( t < tt ) *m++ = *t++;
2109  rdo[1] = WORDDIF(m,rdo);
2110 
2111  *m++ = AC.lUniTrace[0];
2112  *m++ = AC.lUniTrace[1];
2113  *m++ = AC.lUniTrace[2];
2114  *m++ = AC.lUniTrace[3];
2115  *m++ = SNUMBER;
2116  *m++ = 4;
2117  *m++ = 2*norm;
2118  *m++ = -1;
2119 
2120  while ( t < r ) {
2121  if ( *t != GAMMA || t[FUNHEAD] != num ) {
2122  j = t[1];
2123  NCOPY(m,t,j);
2124  }
2125  else t += t[1];
2126  }
2127  rr = term + *term;
2128  while ( t < rr ) *m++ = *t++;
2129 
2130  *termout = WORDDIF(m,termout);
2131  rr = m;
2132  t = termout;
2133  j = *termout;
2134  NCOPY(m,t,j);
2135  AT.WorkPointer = m;
2136  if ( Generator(BHEAD t,level) ) goto ChisCall;
2137 
2138  j = WORDDIF(termout,mat)-1;
2139  t = matpoint;
2140  m = t + j;
2141  AT.WorkPointer = rr;
2142  while ( m > t ) {
2143  i = *--m; *m = *t; *t++ = i;
2144  }
2145 
2146  if ( Generator(BHEAD termout,level) ) goto ChisCall;
2147  AT.WorkPointer = mat;
2148 
2149  goto NextK;
2150  }
2151  t++;
2152  }
2153  }
2154  else t += t[1];
2155  }
2156  s++;
2157  }
2158 NextK:;
2159  }
2160  return(0);
2161 /*
2162  #] Do :
2163 */
2164 ChisCall:
2165  if ( AM.tracebackflag ) {
2166  MLOCK(ErrorMessageLock);
2167  MesCall("Chisholm");
2168  MUNLOCK(ErrorMessageLock);
2169  }
2170  return(-1);
2171 }
2172 
2173 /*
2174  #] Chisholm :
2175  #[ TenVecFind : WORD TenVecFind(term,params)
2176 */
2177 
2178 WORD TenVecFind(PHEAD WORD *term, WORD *params)
2179 {
2180  GETBIDENTITY
2181  WORD *t, *w, *m, *tstop;
2182  WORD i, mode, thevector, thetensor, spectator;
2183  thetensor = params[3];
2184  thevector = params[4];
2185  mode = params[5];
2186  if ( thetensor < 0 ) { /* $-expression */
2187  thetensor = DolToTensor(BHEAD -thetensor);
2188  if ( thetensor < FUNCTION ) {
2189  if ( thevector > 0 ) {
2190  thetensor = DolToTensor(BHEAD thevector);
2191  if ( thetensor < FUNCTION ) {
2192  MLOCK(ErrorMessageLock);
2193  MesPrint("$%s should have been a tensor in module %l"
2194  ,DOLLARNAME(Dollars,params[4]),AC.CModule);
2195  MUNLOCK(ErrorMessageLock);
2196  return(-1);
2197  }
2198  thevector = DolToVector(BHEAD -params[3]);
2199  if ( thevector >= 0 ) {
2200  MLOCK(ErrorMessageLock);
2201  MesPrint("$%s should have been a vector in module %l"
2202  ,DOLLARNAME(Dollars,-params[3]),AC.CModule);
2203  MUNLOCK(ErrorMessageLock);
2204  return(-1);
2205  }
2206  }
2207  else {
2208  MLOCK(ErrorMessageLock);
2209  MesPrint("$%s should have been a tensor in module %l"
2210  ,DOLLARNAME(Dollars,-params[3]),AC.CModule);
2211  MUNLOCK(ErrorMessageLock);
2212  return(-1);
2213  }
2214  }
2215  }
2216  if ( thevector > 0 ) { /* $-expression */
2217  thevector = DolToVector(BHEAD thevector);
2218  if ( thevector >= 0 ) {
2219  MLOCK(ErrorMessageLock);
2220  MesPrint("$%s should have been a vector in module %l"
2221  ,DOLLARNAME(Dollars,params[4]),AC.CModule);
2222  MUNLOCK(ErrorMessageLock);
2223  return(-1);
2224  }
2225  }
2226  if ( ( mode & 1 ) != 0 ) { /* Vector to tensor */
2227  GETSTOP(term,tstop);
2228  t = term + 1;
2229  while ( t < tstop ) {
2230  if ( *t == DOTPRODUCT ) {
2231  i = t[1] - 2; t += 2;
2232  while ( i > 0 ) {
2233  spectator = 0;
2234  if ( t[2] < 0 ) {}
2235  else if ( *t == thevector && t[1] == thevector ) {
2236  if ( ( mode & 2 ) == 0 ) spectator = thevector;
2237  }
2238  else if ( *t == thevector ) spectator = t[1];
2239  else if ( t[1] == thevector ) spectator = *t;
2240  if ( spectator ) {
2241  if ( ( mode & 8 ) == 0 ) goto match;
2242  w = SetElements + Sets[params[6]].first;
2243  m = SetElements + Sets[params[6]].last;
2244  while ( w < m ) {
2245  if ( *w == spectator ) break;
2246  w++;
2247  }
2248  if ( w >= m ) goto match;
2249  }
2250  t += 3;
2251  i -= 3;
2252  }
2253  }
2254  else if ( *t == VECTOR ) {
2255  i = t[1] - 2; t += 2;
2256  while ( i > 0 ) {
2257  if ( *t == thevector ) goto match;
2258  t += 2;
2259  i -= 2;
2260  }
2261  }
2262  else if ( *t == thetensor ) t += t[1];
2263  else if ( *t >= FUNCTION ) {
2264  if ( functions[*t-FUNCTION].spec > 0 ) {
2265  w = t + t[1];
2266  t += FUNHEAD;
2267  while ( t < w ) {
2268  if ( *t == thevector ) goto match;
2269  t++;
2270  }
2271  }
2272  else if ( ( mode & 4 ) != 0 ) {
2273  w = t + t[1];
2274  t += FUNHEAD;
2275  while ( t < w ) {
2276  if ( *t == -VECTOR && t[1] == thevector ) goto match;
2277  else if ( *t > 0 ) t += *t;
2278  else if ( *t <= -FUNCTION ) t++;
2279  else t += 2;
2280  }
2281  }
2282  else t += t[1];
2283  }
2284  else t += t[1];
2285  }
2286  }
2287  else { /* Tensor to Vector */
2288  GETSTOP(term,tstop);
2289  t = term+1;
2290  while ( t < tstop ) {
2291  if ( *t == thetensor ) goto match;
2292  t += t[1];
2293  }
2294  }
2295  return(0);
2296 match:
2297  AT.TMout[0] = 5;
2298  AT.TMout[1] = TENVEC;
2299  AT.TMout[2] = thetensor;
2300  AT.TMout[3] = thevector;
2301  AT.TMout[4] = mode;
2302  if ( ( mode & 8 ) != 0 ) { AT.TMout[0] = 6; AT.TMout[5] = params[6]; }
2303  return(1);
2304 
2305 }
2306 
2307 /*
2308  #] TenVecFind :
2309  #[ TenVec : WORD TenVec(term,params,num,level)
2310 */
2311 
2312 WORD TenVec(PHEAD WORD *term, WORD *params, WORD num, WORD level)
2313 {
2314  GETBIDENTITY
2315  WORD *t, *m, *w, *termout, *tstop, *outlist, *ou, *ww, *mm;
2316  WORD i, j, k, x, mode, thevector, thetensor, DumNow, spectator;
2317  DUMMYUSE(num);
2318  thetensor = params[2];
2319  thevector = params[3];
2320  mode = params[4];
2321  termout = AT.WorkPointer;
2322  DumNow = AR.CurDum = DetCurDum(BHEAD term);
2323  if ( ( mode & 1 ) != 0 ) { /* Vector to tensor */
2324  AT.WorkPointer += *term;
2325  ou = outlist = AT.WorkPointer;
2326  GETSTOP(term,tstop);
2327  t = term + 1;
2328  m = termout + 1;
2329  while ( t < tstop ) {
2330  if ( *t == DOTPRODUCT ) {
2331  i = t[1] - 2;
2332  w = m;
2333  *m++ = *t++; *m++ = *t++;
2334  while ( i > 0 ) {
2335  spectator = 0;
2336  if ( t[2] < 0 ) {
2337  *m++ = *t++; *m++ = *t++; *m++ = *t++;
2338  }
2339  else if ( *t == thevector && t[1] == thevector ) {
2340  if ( ( mode & 2 ) == 0 ) spectator = thevector;
2341  else {
2342  *m++ = *t++; *m++ = *t++; *m++ = *t++;
2343  }
2344  }
2345  else if ( *t == thevector ) spectator = t[1];
2346  else if ( t[1] == thevector ) spectator = *t;
2347  else {
2348  *m++ = *t++; *m++ = *t++; *m++ = *t++;
2349  }
2350  if ( spectator ) {
2351  if ( ( mode & 8 ) == 0 ) goto noveto;
2352  ww = SetElements + Sets[params[5]].first;
2353  mm = SetElements + Sets[params[5]].last;
2354  while ( ww < mm ) {
2355  if ( *ww == spectator ) break;
2356  ww++;
2357  }
2358  if ( ww < mm ) {
2359  *m++ = *t++; *m++ = *t++; *m++ = *t++;
2360  }
2361  else {
2362 noveto: if ( spectator == thevector ) {
2363  for ( j = 0; j < t[2]; j++ ) {
2364  *ou++ = ++AR.CurDum;
2365  *ou++ = AR.CurDum;
2366  }
2367  t += 3;
2368  }
2369  else {
2370  for ( j = 0; j < t[2]; j++ ) *ou++ = spectator;
2371  t += 3;
2372  }}
2373  }
2374  i -= 3;
2375  }
2376  w[1] = WORDDIF(m,w);
2377  if ( w[1] == 2 ) m = w;
2378  }
2379  else if ( *t == VECTOR ) {
2380  i = t[1] - 2; w = m;
2381  *m++ = *t++; *m++ = *t++;
2382  while ( i > 0 ) {
2383  if ( *t == thevector ) {
2384  *ou++ = t[1];
2385  t += 2;
2386  }
2387  else { *m++ = *t++; *m++ = *t++; }
2388  i -= 2;
2389  }
2390  w[1] = WORDDIF(m,w);
2391  if ( w[1] == 2 ) m = w;
2392  }
2393  else if ( *t == thetensor ) {
2394  i = t[1] - FUNHEAD;
2395  t += FUNHEAD;
2396  NCOPY(ou,t,i);
2397  }
2398  else if ( *t >= FUNCTION ) {
2399  if ( functions[*t-FUNCTION].spec > 0 ) {
2400  w = t + t[1];
2401  i = FUNHEAD;
2402  NCOPY(m,t,i);
2403  while ( t < w ) {
2404  if ( *t == thevector ) {
2405  *m++ = ++AR.CurDum;
2406  *ou++ = AR.CurDum;
2407  t++;
2408  }
2409  else *m++ = *t++;
2410  }
2411  }
2412  else if ( ( mode & 4 ) != 0 ) {
2413  w = t + t[1];
2414  i = FUNHEAD;
2415  NCOPY(m,t,i);
2416  while ( t < w ) {
2417  if ( *t == -VECTOR && t[1] == thevector ) {
2418  *m++ = -INDEX;
2419  *m++ = ++AR.CurDum;
2420  *ou++ = AR.CurDum;
2421  t += 2;
2422  }
2423  else if ( *t > 0 ) {
2424  i = *t;
2425  NCOPY(m,t,i);
2426  }
2427  else if ( *t <= -FUNCTION ) *m++ = *t++;
2428  else { *m++ = *t++; *m++ = *t++; }
2429  }
2430  }
2431  else goto docopy;
2432  }
2433  else {
2434 docopy:
2435  i = t[1];
2436  NCOPY(m,t,i);
2437  }
2438  }
2439  i = WORDDIF(ou,outlist);
2440  if ( i > 0 ) {
2441  for ( j = 1; j < i; j++ ) {
2442  if ( outlist[j-1] > outlist[j] ) {
2443  x = outlist[j-1]; outlist[j-1] = outlist[j]; outlist[j] = x;
2444  for ( k = j-1; k > 0; k-- ) {
2445  if ( outlist[k-1] <= outlist[k] ) break;
2446  x = outlist[k-1]; outlist[k-1] = outlist[k]; outlist[k] = x;
2447  }
2448  }
2449  }
2450 
2451  *m++ = thetensor;
2452  *m++ = FUNHEAD + i;
2453  *m++ = DIRTYSYMFLAG;
2454  FILLFUN3(m)
2455  ou = outlist;
2456  NCOPY(m,ou,i);
2457  }
2458  w = term + *term;
2459  while ( t < w ) *m++ = *t++;
2460  }
2461  else { /* Tensor to Vector */
2462  GETSTOP(term,tstop);
2463  t = term+1;
2464  m = termout+1;
2465  while ( t < tstop ) {
2466  if ( *t != thetensor ) {
2467  i = t[1];
2468  NCOPY(m,t,i);
2469  }
2470  else {
2471  i = t[1] - FUNHEAD;
2472  t += FUNHEAD;
2473  if ( i > 0 ) {
2474  w = m; m += 2;
2475  while ( --i >= 0 ) {
2476  *m++ = thevector;
2477  *m++ = *t++;
2478  }
2479  *w = DELTA;
2480  w[1] = WORDDIF(m,w);
2481  }
2482  }
2483  }
2484  w = term + *term;
2485  while ( t < w ) *m++ = *t++;
2486  }
2487  *termout = WORDDIF(m,termout);
2488  AT.WorkPointer = m;
2489  *AT.TMout = 0;
2490  if ( Generator(BHEAD termout,level) ) goto fromTenVec;
2491  AR.CurDum = DumNow;
2492  AT.WorkPointer = termout;
2493  return(0);
2494 fromTenVec:
2495  if ( AM.tracebackflag ) {
2496  MLOCK(ErrorMessageLock);
2497  MesCall("TenVec");
2498  MUNLOCK(ErrorMessageLock);
2499  }
2500  return(-1);
2501 }
2502 
2503 /*
2504  #] TenVec :
2505  #] Operations :
2506 */
WORD ** lhs
Definition: structs.h:942
Definition: structs.h:938
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:3101