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