1 /** @file opera.c
2 *
3 * Contains the 'operations'
4 * These are the trace routines, the contractions of the Levi-Civita tensors
5 * and the tensor to vector/vector to tensor routines.
6 * The trace and contraction routines are done in a special way
7 * (see commentary with the FIXEDGLOBALS struct)
8 */
9 /* #[ License : */
10 /*
11 * Copyright (C) 1984-2017 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
EpfFind(PHEAD WORD * term,WORD * params)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
EpfCon(PHEAD WORD * term,WORD * params,WORD num,WORD level)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
EpfGen(WORD number,WORD * inlist,WORD * kron,WORD * perm,WORD sgn)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
Trick(WORD * in,TRACES * t)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
Trace4no(WORD number,WORD * kron,TRACES * t)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
Trace4(PHEAD WORD * term,WORD * params,WORD num,WORD level)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
Trace4Gen(PHEAD TRACES * t,WORD number)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
TraceNno(WORD number,WORD * kron,TRACES * t)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
TraceN(PHEAD WORD * term,WORD * params,WORD num,WORD level)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
TraceNgen(PHEAD TRACES * t,WORD number)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
Traces(PHEAD WORD * term,WORD * params,WORD num,WORD level)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
TraceFind(PHEAD WORD * term,WORD * params)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
Chisholm(PHEAD WORD * term,WORD level)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
TenVecFind(PHEAD WORD * term,WORD * params)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
TenVec(PHEAD WORD * term,WORD * params,WORD num,WORD level)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 */
2507