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