FORM v5.0.1-33-gdf7fc94
lus.c
Go to the documentation of this file.
1
7/* #[ License : */
8/*
9 * Copyright (C) 1984-2026 J.A.M. Vermaseren
10 * When using this file you are requested to refer to the publication
11 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
12 * This is considered a matter of courtesy as the development was paid
13 * for by FOM the Dutch physics granting agency and we would like to
14 * be able to track its scientific use to convince FOM of its value
15 * for the community.
16 *
17 * This file is part of FORM.
18 *
19 * FORM is free software: you can redistribute it and/or modify it under the
20 * terms of the GNU General Public License as published by the Free Software
21 * Foundation, either version 3 of the License, or (at your option) any later
22 * version.
23 *
24 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
25 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
26 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
27 * details.
28 *
29 * You should have received a copy of the GNU General Public License along
30 * with FORM. If not, see <http://www.gnu.org/licenses/>.
31 */
32/* #] License : */
33/*
34 #[ Includes : lus.c
35*/
36
37#include "form3.h"
38
39/*
40 #] Includes :
41 #[ Lus :
42
43 Routine to find loops.
44 Mode: 0: Just tell whether there is such a loop.
45 1: Take out the functions and replace by outfun with the
46 remaining arguments of the loop function
47 > AM.OffsetIndex: This index must be included in the loop.
48 < -AM.OffsetIndex: This index must be included in the loop. Replace.
49 Return value: 0: no loop. 1: there is/was such a loop.
50 funnum: the function(s) in which we look for a loop.
51 numargs: the number of arguments admissible in the function.
52 outfun: the output function in case of a substitution.
53 loopsize: the size of the loop we are looking for.
54 if < 0 we look for all loops.
55*/
56
57int Lus(WORD *term, WORD funnum, WORD loopsize, WORD numargs, WORD outfun, WORD mode)
58{
59 GETIDENTITY
60 WORD *w, *t, *tt, *m, *r, **loc, *tstop, minloopsize;
61 int nfun, i, j, jj, k, n, sign = 0, action = 0, L, ten, ten2, totnum,
62 sign2, *alist, *wi, mini, maxi, medi = 0;
63 if ( numargs <= 1 ) return(0);
64 GETSTOP(term,tstop);
65/*
66 First count the number of functions with the proper number of arguments.
67*/
68 t = term+1; nfun = 0;
69 if ( ( ten = functions[funnum-FUNCTION].spec ) >= TENSORFUNCTION ) {
70 while ( t < tstop ) {
71 if ( *t == funnum && t[1] == FUNHEAD+numargs ) { nfun++; }
72 t += t[1];
73 }
74 }
75 else {
76 while ( t < tstop ) {
77 if ( *t == funnum ) {
78 i = 0; m = t+FUNHEAD; t += t[1];
79 while ( m < t ) { i++; NEXTARG(m) }
80 if ( i == numargs ) nfun++;
81 }
82 else t += t[1];
83 }
84 }
85 if ( loopsize < 0 ) minloopsize = 2;
86 else minloopsize = loopsize;
87 if ( funnum < minloopsize ) return(0); /* quick abort */
88 if ( ((functions[funnum-FUNCTION].symmetric) & ~REVERSEORDER) == ANTISYMMETRIC ) sign = 1;
89 if ( mode == 1 || mode < 0 ) {
90 ten2 = functions[outfun-FUNCTION].spec >= TENSORFUNCTION;
91 }
92 else ten2 = -1;
93/*
94 Allocations:
95*/
96 if ( AN.numflocs < funnum ) {
97 if ( AN.funlocs ) M_free(AN.funlocs,"Lus: AN.funlocs");
98 AN.numflocs = funnum;
99 AN.funlocs = (WORD **)Malloc1(sizeof(WORD *)*AN.numflocs,"Lus: AN.funlocs");
100 }
101 if ( AN.numfargs < funnum*numargs ) {
102 if ( AN.funargs ) M_free(AN.funargs,"Lus: AN.funargs");
103 AN.numfargs = funnum*numargs;
104 AN.funargs = (int *)Malloc1(sizeof(int *)*AN.numfargs,"Lus: AN.funargs");
105 }
106/*
107 Make a list of relevant indices
108*/
109 alist = AN.funargs; loc = AN.funlocs;
110 t = term+1;
111 if ( ten >= TENSORFUNCTION ) {
112 while ( t < tstop ) {
113 if ( *t == funnum && t[1] == FUNHEAD+numargs ) {
114 *loc++ = t;
115 t += FUNHEAD;
116 j = i = numargs; while ( --i >= 0 ) {
117 if ( *t >= AM.OffsetIndex &&
118 ( *t >= AM.OffsetIndex+WILDOFFSET ||
119 indices[*t-AM.OffsetIndex].dimension != 0 ) ) {
120 *alist++ = *t++; j--;
121 }
122 else t++;
123 }
124 while ( --j >= 0 ) *alist++ = -1;
125 }
126 else t += t[1];
127 }
128 }
129 else {
130 nfun = 0;
131 while ( t < tstop ) {
132 if ( *t == funnum ) {
133 w = t;
134 i = 0; m = t+FUNHEAD; t += t[1];
135 while ( m < t ) { i++; NEXTARG(m) }
136 if ( i == numargs ) {
137 m = w + FUNHEAD;
138 while ( m < t ) {
139 if ( *m == -INDEX && m[1] >= AM.OffsetIndex &&
140 ( m[1] >= AM.OffsetIndex+WILDOFFSET ||
141 indices[m[1]-AM.OffsetIndex].dimension != 0 ) ) {
142 *alist++ = m[1]; m += 2; i--;
143 }
144 else if ( ten2 >= TENSORFUNCTION && *m != -INDEX
145 && *m != -VECTOR && *m != -MINVECTOR &&
146 ( *m != -SNUMBER || *m < 0 || *m >= AM.OffsetIndex ) ) {
147 i = numargs; break;
148 }
149 else { NEXTARG(m) }
150 }
151 if ( i < numargs ) {
152 *loc++ = w;
153 nfun++;
154 while ( --i >= 0 ) *alist++ = -1;
155 }
156 }
157 }
158 else t += t[1];
159 }
160 if ( nfun < minloopsize ) return(0);
161 }
162/*
163 We have now nfun objects. Not all indices may be usable though.
164 If the list is not long, we use a quadratic algorithm to remove
165 indices and vertices that cannot be used. If it becomes large we
166 sort the list of available indices (and their multiplicity) and
167 work with binary searches.
168*/
169 alist = AN.funargs; totnum = numargs*nfun;
170 if ( nfun > 7 ) {
171 if ( AN.funisize < totnum ) {
172 if ( AN.funinds ) M_free(AN.funinds,"AN.funinds");
173 AN.funisize = (totnum*3)/2;
174 AN.funinds = (int *)Malloc1(AN.funisize*2*sizeof(int),"AN.funinds");
175 }
176 i = totnum; n = 0; wi = AN.funinds;
177 while ( --i >= 0 ) {
178 if ( *alist >= 0 ) { n++; *wi++ = *alist; *wi++ = 1; }
179 alist++;
180 }
181 n = SortTheList(AN.funinds,n);
182 do {
183 action = 0;
184 for ( i = 0; i < nfun; i++ ) {
185 alist = AN.funargs + i*numargs;
186 jj = numargs;
187 for ( j = 0; j < jj; j++ ) {
188 if ( alist[j] < 0 ) break;
189 mini = 0; maxi = n-1;
190 while ( mini <= maxi ) {
191 medi = (mini + maxi) / 2; k = AN.funinds[2*medi];
192 if ( alist[j] > k ) mini = medi + 1;
193 else if ( alist[j] < k ) maxi = medi - 1;
194 else break;
195 }
196 if ( AN.funinds[2*medi+1] <= 1 ) {
197 (AN.funinds[2*medi+1])--;
198 jj--; k = j; while ( k < jj ) { alist[k] = alist[k+1]; k++; }
199 alist[jj] = -1; j--;
200 }
201 }
202 if ( jj < 2 ) {
203 if ( jj == 1 ) {
204 mini = 0; maxi = n-1;
205 while ( mini <= maxi ) {
206 medi = (mini + maxi) / 2; k = AN.funinds[2*medi];
207 if ( alist[0] > k ) mini = medi + 1;
208 else if ( alist[0] < k ) maxi = medi - 1;
209 else break;
210 }
211 (AN.funinds[2*medi+1])--;
212 if ( AN.funinds[2*medi+1] == 1 ) action++;
213 }
214 nfun--; totnum -= numargs; AN.funlocs[i] = AN.funlocs[nfun];
215 wi = AN.funargs + nfun*numargs;
216 for ( j = 0; j < numargs; j++ ) alist[j] = *wi++;
217 i--;
218 }
219 }
220 } while ( action );
221 }
222 else {
223 for ( i = 0; i < totnum; i++ ) {
224 if ( alist[i] == -1 ) continue;
225 for ( j = 0; j < totnum; j++ ) {
226 if ( alist[j] == alist[i] && j != i ) break;
227 }
228 if ( j >= totnum ) alist[i] = -1;
229 }
230 do {
231 action = 0;
232 for ( i = 0; i < nfun; i++ ) {
233 alist = AN.funargs + i*numargs;
234 n = numargs;
235 for ( k = 0; k < n; k++ ) {
236 if ( alist[k] < 0 ) { alist[k--] = alist[--n]; alist[n] = -1; }
237 }
238 if ( n <= 1 ) {
239 if ( n == 1 ) { j = alist[0]; }
240 else j = -1;
241 nfun--; totnum -= numargs; AN.funlocs[i] = AN.funlocs[nfun];
242 wi = AN.funargs + nfun * numargs;
243 for ( k = 0; k < numargs; k++ ) alist[k] = wi[k];
244 i--;
245 if ( j >= 0 ) {
246 for ( k = 0, jj = 0, wi = AN.funargs; k < totnum; k++, wi++ ) {
247 if ( *wi == j ) { jj++; if ( jj > 1 ) break; }
248 }
249 if ( jj <= 1 ) {
250 for ( k = 0, wi = AN.funargs; k < totnum; k++, wi++ ) {
251 if ( *wi == j ) { *wi = -1; action = 1; }
252 }
253 }
254 }
255 }
256 }
257 } while ( action );
258 }
259 if ( nfun < minloopsize ) return(0);
260/*
261 Now we have nfun objects, each with at least 2 indices, each of which
262 occurs at least twice in our list. There will be a loop!
263*/
264 if ( mode != 0 && mode != 1 ) {
265 if ( mode > 0 ) AN.tohunt = mode - 5;
266 else AN.tohunt = -mode - 5;
267 AN.nargs = numargs; AN.numoffuns = nfun;
268 i = 0;
269 if ( loopsize < 0 ) {
270 if ( loopsize == -1 ) k = nfun;
271 else { k = -loopsize-1; if ( k > nfun ) k = nfun; }
272 for ( L = 2; L <= k; L++ ) {
273 if ( FindLus(0,L,AN.tohunt) ) goto Success;
274 }
275 }
276 else if ( FindLus(0,loopsize,AN.tohunt) ) { L = loopsize; goto Success; }
277 }
278 else {
279 AN.nargs = numargs; AN.numoffuns = nfun;
280 if ( loopsize < 0 ) {
281 jj = 2; k = nfun;
282 if ( loopsize < -1 ) { k = -loopsize-1; if ( k > nfun ) k = nfun; }
283 }
284 else { jj = k = loopsize; }
285 for ( L = jj; L <= k; L++ ) {
286 for ( i = 0; i <= nfun-L; i++ ) {
287 alist = AN.funargs + i * numargs;
288 for ( jj = 0; jj < numargs; jj++ ) {
289 if ( alist[jj] < 0 ) continue;
290 AN.tohunt = alist[jj];
291 for ( j = jj+1; j < numargs; j++ ) {
292 if ( alist[j] < 0 ) continue;
293 if ( FindLus(i+1,L-1,alist[j]) ) {
294 alist[0] = alist[jj];
295 alist[1] = alist[j];
296 goto Success;
297 }
298 }
299 }
300 }
301 }
302 }
303 return(0);
304Success:;
305 if ( mode == 0 || mode > 1 ) return(1);
306/*
307 Now we have to make the replacement and fix the potential sign
308*/
309 sign2 = 1;
310 wi = AN.funargs + i*numargs; loc = AN.funlocs + i;
311 for ( k = 0; k < L; k++ ) *(loc[k]) = -1;
312 if ( AT.WorkPointer < term + *term ) AT.WorkPointer = term + *term;
313 w = AT.WorkPointer + 1;
314 m = t = term + 1;
315 while ( t < tstop ) {
316 if ( *t == -1 ) break;
317 t += t[1];
318 }
319 while ( m < t ) *w++ = *m++;
320 r = w;
321 *w++ = outfun;
322 w++;
323 *w++ = DIRTYFLAG;
324 FILLFUN3(w)
325 if ( functions[outfun-FUNCTION].spec >= TENSORFUNCTION ) {
326 if ( ten >= TENSORFUNCTION ) {
327 for ( i = 0; i < L; i++ ) {
328 alist = wi + i*numargs;
329 m = loc[i] + FUNHEAD;
330 for ( k = 0; k < numargs; k++ ) {
331 if ( m[k] == alist[0] ) {
332 if ( k != 0 ) {
333 jj = m[k]; m[k] = m[0]; m[0] = jj;
334 sign = -sign;
335 }
336 break;
337 }
338 }
339 for ( k = 1; k < numargs; k++ ) {
340 if ( m[k] == alist[1] ) {
341 if ( k != 1 ) {
342 jj = m[k]; m[k] = m[1]; m[1] = jj;
343 sign = -sign;
344 }
345 break;
346 }
347 }
348 m += 2;
349 for ( k = 2; k < numargs; k++ ) *w++ = *m++;
350 }
351 }
352 else {
353 WORD *t1, *t2, *t3;
354 for ( i = 0; i < L; i++ ) {
355 alist = wi + i*numargs;
356 tt = loc[i];
357 m = tt + FUNHEAD;
358 for ( k = 0; k < numargs; k++ ) {
359 if ( *m == -INDEX && m[1] == alist[0] ) {
360 if ( k != 0 ) {
361 if ( ( k & 1 ) != 0 ) sign = -sign;
362/*
363 now move to position 0
364*/
365 t2 = m+2; t1 = m; t3 = tt+FUNHEAD;
366 while ( t1 > t3 ) { *--t2 = *--t1; }
367 t3[0] = -INDEX; t3[1] = alist[0];
368 }
369 break;
370 }
371 NEXTARG(m)
372 }
373 m = tt + FUNHEAD + 2;
374 for ( k = 1; k < numargs; k++ ) {
375 if ( *m == -INDEX && m[1] == alist[1] ) {
376 if ( k != 1 ) {
377 if ( ( k & 1 ) == 0 ) sign = -sign;
378/*
379 now move to position 1
380*/
381 t2 = m+2; t1 = m; t3 = tt+FUNHEAD+2;
382 while ( t1 > t3 ) { *--t2 = *--t1; }
383 t3[0] = -INDEX; t3[1] = alist[1];
384 }
385 break;
386 }
387 NEXTARG(m)
388 }
389/*
390 now copy the remaining arguments to w
391 keep in mind that the output function is a tensor!
392*/
393 t1 = tt + FUNHEAD + 4;
394 t2 = tt + tt[1];
395 while ( t1 < t2 ) {
396 if ( *t1 == -INDEX || *t1 == -VECTOR ) {
397 *w++ = t1[1]; t1 += 2;
398 }
399 else if ( *t1 == -MINVECTOR ) {
400 *w++ = t1[1]; t1 += 2; sign2 = -sign2;
401 }
402 else if ( ( *t1 == -SNUMBER ) && ( t1[1] >= 0 ) && ( t1[1] < AM.OffsetIndex ) ) {
403 *w++ = t1[1]; t1 += 2; sign2 = -sign2;
404 }
405 else {
406 MLOCK(ErrorMessageLock);
407 MesPrint("Illegal attempt to use a non-index-like argument in a tensor in ReplaceLoop statement");
408 MUNLOCK(ErrorMessageLock);
409 Terminate(-1);
410 }
411 }
412 }
413 }
414 }
415 else {
416 if ( ten >= TENSORFUNCTION ) {
417 for ( i = 0; i < L; i++ ) {
418 alist = wi + i*numargs;
419 m = loc[i] + FUNHEAD;
420 for ( k = 0; k < numargs; k++ ) {
421 if ( m[k] == alist[0] ) {
422 if ( k != 0 ) {
423 jj = m[k]; m[k] = m[0]; m[0] = jj;
424 sign = -sign;
425 break;
426 }
427 }
428 }
429 for ( k = 1; k < numargs; k++ ) {
430 if ( m[k] == alist[1] ) {
431 if ( k != 1 ) {
432 jj = m[k]; m[k] = m[1]; m[1] = jj;
433 sign = -sign;
434 break;
435 }
436 }
437 }
438 m += 2;
439 for ( k = 2; k < numargs; k++ ) {
440 if ( *m >= AM.OffsetIndex ) { *w++ = -INDEX; }
441 else if ( *m < 0 ) { *w++ = -VECTOR; }
442 else { *w = -SNUMBER; }
443 *w++ = *m++;
444 }
445 }
446 }
447 else {
448 WORD *t1, *t2, *t3;
449 for ( i = 0; i < L; i++ ) {
450 alist = wi + i*numargs;
451 tt = loc[i];
452 m = tt + FUNHEAD;
453 for ( k = 0; k < numargs; k++ ) {
454 if ( *m == -INDEX && m[1] == alist[0] ) {
455 if ( k != 0 ) {
456 if ( ( k & 1 ) != 0 ) sign = -sign;
457/*
458 now move to position 0
459*/
460 t2 = m+2; t1 = m; t3 = tt+FUNHEAD;
461 while ( t1 > t3 ) { *--t2 = *--t1; }
462 t3[0] = -INDEX; t3[1] = alist[0];
463 }
464 break;
465 }
466 NEXTARG(m)
467 }
468 m = tt + FUNHEAD + 2;
469 for ( k = 1; k < numargs; k++ ) {
470 if ( *m == -INDEX && m[1] == alist[1] ) {
471 if ( k != 1 ) {
472 if ( ( k & 1 ) == 0 ) sign = -sign;
473/*
474 now move to position 1
475*/
476 t2 = m+2; t1 = m; t3 = tt+FUNHEAD+2;
477 while ( t1 > t3 ) { *--t2 = *--t1; }
478 t3[0] = -INDEX; t3[1] = alist[1];
479 }
480 break;
481 }
482 NEXTARG(m)
483 }
484/*
485 now copy the remaining arguments to w
486*/
487 t1 = tt + FUNHEAD + 4;
488 t2 = tt + tt[1];
489 while ( t1 < t2 ) *w++ = *t1++;
490 }
491 }
492 }
493 r[1] = w-r;
494 while ( t < tstop ) {
495 if ( *t == -1 ) { t += t[1]; continue; }
496 i = t[1];
497 NCOPY(w,t,i)
498 }
499 tstop = term + *term;
500 while ( t < tstop ) *w++ = *t++;
501 if ( sign < 0 ) w[-1] = -w[-1];
502 i = w - AT.WorkPointer;
503 *AT.WorkPointer = i;
504 t = term; w = AT.WorkPointer;
505 NCOPY(t,w,i)
506 *AN.RepPoint = 1; /* For Repeat */
507 return(1);
508}
509
510/*
511 #] Lus :
512 #[ FindLus :
513*/
514
515int FindLus(int from, int level, int openindex)
516{
517 GETIDENTITY
518 int i, j, k, jj, *alist, *blist, *w, *m, partner;
519 WORD **loc = AN.funlocs, *wor;
520 if ( level == 1 ) {
521 for ( i = from; i < AN.numoffuns; i++ ) {
522 alist = AN.funargs + i*AN.nargs;
523 for ( j = 0; j < AN.nargs; j++ ) {
524 if ( alist[j] == openindex ) {
525 for ( k = 0; k < AN.nargs; k++ ) {
526 if ( k == j ) continue;
527 if ( alist[k] == AN.tohunt ) {
528 loc[from] = loc[i];
529 alist = AN.funargs + from*AN.nargs;
530 alist[0] = openindex; alist[1] = AN.tohunt;
531 return(1);
532 }
533 }
534 }
535 }
536 }
537 }
538 else {
539 for ( i = from; i < AN.numoffuns; i++ ) {
540 alist = AN.funargs + i*AN.nargs;
541 for ( j = 0; j < AN.nargs; j++ ) {
542 if ( alist[j] == openindex ) {
543 if ( from != i ) {
544 wor = loc[i]; loc[i] = loc[from]; loc[from] = wor;
545 blist = w = AN.funargs + from*AN.nargs;
546 m = alist;
547 k = AN.nargs;
548 while ( --k >= 0 ) { jj = *m; *m++ = *w; *w++ = jj; }
549 }
550 else blist = alist;
551 for ( k = 0; k < AN.nargs; k++ ) {
552 if ( k == j || blist[k] < 0 ) continue;
553 partner = blist[k];
554 if ( FindLus(from+1,level-1,partner) ) {
555 blist[0] = openindex; blist[1] = partner;
556 return(1);
557 }
558 }
559 if ( from != i ) {
560 wor = loc[i]; loc[i] = loc[from]; loc[from] = wor;
561 w = AN.funargs + from*AN.nargs;
562 m = alist;
563 k = AN.nargs;
564 while ( --k >= 0 ) { jj = *m; *m++ = *w; *w++ = jj; }
565 }
566 }
567 }
568 }
569 }
570 return(0);
571}
572
573/*
574 #] FindLus :
575 #[ SortTheList :
576*/
577
578int SortTheList(int *slist, int num)
579{
580 GETIDENTITY
581 int i, nleft, nright, *t1, *t2, *t3, *rlist;
582 if ( num <= 2 ) {
583 if ( num <= 1 ) return(num);
584 if ( slist[0] < slist[2] ) return(2);
585 if ( slist[0] > slist[2] ) {
586 i = slist[0]; slist[0] = slist[2]; slist[2] = i;
587 i = slist[1]; slist[1] = slist[3]; slist[3] = i;
588 return(2);
589 }
590 slist[1] += slist[3];
591 return(1);
592 }
593 else {
594 nleft = num/2; rlist = slist + 2*nleft;
595 nright = SortTheList(rlist,num-nleft);
596 nleft = SortTheList(slist,nleft);
597 if ( AN.tlistsize < nleft ) {
598 if ( AN.tlistbuf ) M_free(AN.tlistbuf,"AN.tlistbuf");
599 AN.tlistsize = (nleft*3)/2;
600 AN.tlistbuf = (int *)Malloc1(AN.tlistsize*2*sizeof(int),"AN.tlistbuf");
601 }
602 i = nleft; t1 = slist; t2 = AN.tlistbuf;
603 while ( --i >= 0 ) { *t2++ = *t1++; *t2++ = *t1++; }
604 i = nleft+nright; t1 = AN.tlistbuf; t2 = rlist; t3 = slist;
605 while ( nleft > 0 && nright > 0 ) {
606 if ( *t1 < *t2 ) {
607 *t3++ = *t1++; *t3++ = *t1++; nleft--;
608 }
609 else if ( *t1 > *t2 ) {
610 *t3++ = *t2++; *t3++ = *t2++; nright--;
611 }
612 else {
613 *t3++ = *t1++; t2++; *t3++ = (*t1++) + (*t2++); i--;
614 nleft--; nright--;
615 }
616 }
617 while ( --nleft >= 0 ) { *t3++ = *t1++; *t3++ = *t1++; }
618 while ( --nright >= 0 ) { *t3++ = *t2++; *t3++ = *t2++; }
619 return(i);
620 }
621}
622
623/*
624 #] SortTheList :
625 #[ AllLoops :
626
627 Routine finds all possible loops that can be made in the
628 arguments of the symmetric commuting vertex function v and creates
629 for each loop a new term which is the original term times the loop.
630 The occurrences of v can have different numbers of arguments.
631 Each argument that occurs twice (and in different instances of v)
632 in total will be considered.
633
634 The input parameters are in C->lhs[level] and are
635 C->lhs[level][0] TYPEALLLOOPS
636 C->lhs[level][1] 6
637 C->lhs[level][2] Number of function v.
638 C->lhs[level][3] Number of the output function loop.
639 C->lhs[level][4] option1: the type of argument.
640 C->lhs[level][5] option2: 0,1: what to do if no loop.
641 Eligible arguments must be of the type SYMBOL, VECTOR, INDEX or SNUMBER.
642 They must all be of the same type, indicated by option1.
643 Extra restriction: In a loop, each v can be visited at most once.
644 If there is no loop at all, option2 determines whether the term
645 remains unmodified, or is replaced by zero.
646 Function v can be either a regular function or a tensor.
647 To facilitate this we copy the relevant arguments into the workspace.
648*/
649
650int AllLoops(PHEAD WORD *term,WORD level)
651{
652 CBUF *C = cbuf+AM.rbufnum;
653 WORD vcode = C->lhs[level][2]; /* The input function */
654 WORD option1 = C->lhs[level][4]; /* type of argument */
655 WORD option2 = C->lhs[level][5]; /* what to do when no loop */
656 WORD *tstop, *t, *tend, *tstart, *tfrom;
657 WORD i, j, jj;
658 WORD *arglist, nargs, *loop, nloop;
659 WORD *oldworkpointer = AT.WorkPointer;
660 LONG oldpworkpointer = AT.pWorkPointer;
661 LONG numgenerated = 0, vert;
662 WORD *a, *a1, *a2, *a3, *v, *vv, nvert, *to, *from, *tos, action;
663/*
664 Search for the first occurrence of vcode.
665*/
666 tstop = term+*term; tstop -= ABS(tstop[-1]);
667 t = term + 1;
668 while ( t < tstop && *t != vcode ) t += t[1];
669 if ( t == tstop ) {
670 if ( option2 == 0 ) return(0);
671 else return(Generator(BHEAD term,level));
672 }
673 tstart = t;
674 nvert = 0;
675 do {
676 t += t[1]; nvert++;
677 } while ( t < tstop && *t == vcode );
678 tend = t;
679/*
680 Make room for 2*nvert pointers
681*/
682 WantAddPointers(2*nvert);
683 vert = AT.pWorkPointer;
684 AT.pWorkPointer += 2*nvert;
685 nvert = 0;
686/*
687 Next we copy these functions into the workspace, but only the arguments
688 that agree with option1. Because of the difference between tensors and
689 regular functions in the copy we strip the type of the argument.
690*/
691 to = AT.WorkPointer;
692 from = tstart;
693 if ( functions[vcode-FUNCTION].spec == TENSORFUNCTION ) {
694 from = tstart;
695 while ( from < tend ) {
696 tos = to++;
697 tfrom = from+from[1];
698 from += FUNHEAD;
699 while ( from < tfrom ) {
700 if ( option1 == -INDEX && *from >= 0 ) {
701 *to++ = *from++;
702 }
703 else if ( option1 == -VECTOR && *from < 0 ) {
704 *to++ = *from++ - AM.OffsetVector;
705 }
706 else {
707 from++;
708 }
709 }
710 *tos = to-from;
711 if ( *tos < 3 ) to = tos;
712 else AT.pWorkSpace[vert+nvert++] = tos;
713 }
714 }
715 else if ( functions[vcode-FUNCTION].spec == 0 ) {
716 from = tstart;
717 while ( from < tend ) {
718 tos = to++;
719 tfrom = from + from[1];
720 from += FUNHEAD;
721 while ( from < tfrom ) {
722 if ( option1 == -VECTOR
723 && ( from[0] == -INDEX || from[0] == -VECTOR )
724 && from[1] < 0 ) {
725 from++;
726 *to++ = *from++ - AM.OffsetVector;
727 }
728 else if ( option1 == *from ) {
729 from++; *to++ = *from++;
730 }
731 else {
732 NEXTARG(from);
733 }
734 }
735 *tos = to-tos;
736 if ( *tos < 3 ) to = tos;
737 else AT.pWorkSpace[vert+nvert++] = tos;
738 }
739 }
740 else {
741 AT.WorkPointer = oldworkpointer;
742 AT.pWorkPointer = oldpworkpointer;
743 if ( option2 == 0 ) return(0);
744 return(Generator(BHEAD term,level));
745 }
746 AT.WorkPointer = to;
747/*
748 Now make a list of all relevant arguments.
749*/
750 a = arglist = AT.WorkPointer;
751 for ( i = 0; i < nvert; i++ ) {
752 v = AT.pWorkSpace[vert+i];
753 j = *v-1; v++;
754 NCOPY(a,v,j);
755 }
756 nargs = a-arglist;
757 AT.WorkPointer = a;
758/*
759 Now sort the list. Bubble type sort.
760*/
761 a1 = arglist; a2 = a1+1;
762 while ( a2 < a ) {
763 a3 = a2;
764 while ( a3 > a1 && a3[-1] > a3[0] ) {
765 EXCH(*a3,a3[-1]);
766 a3--;
767 }
768 a2++;
769 }
770/*
771 Now remove the elements from the list that do not occur exactly twice.
772*/
773 a1 = arglist; a2 = a1; a3 = a1+nargs;
774 while ( a2 < a3 ) {
775 if ( a2+1 == a3 ) { break; }
776 else if ( a2+2 == a3 ) {
777 if ( a2[0] == a2[1] ) { *a1++ = a2[0]; }
778 break;
779 }
780 else {
781 if ( a2[0] != a2[1] ) { a2++; }
782 else if ( a2[0] != a2[2] ) {
783 *a1++ = a2[0]; a2 += 2;
784 }
785 else {
786 a2++; while ( a2 < a3 && a2[-1] == a2[0] ) a2++;
787 }
788 }
789 }
790 nargs = a1-arglist;
791/*
792 Now we need to redo the list of vertices and remove the elements
793 that are not in our arglist.
794*/
795 do {
796 action = 0;
797 for ( i = 0; i < nvert; i++ ) {
798 vv = v = AT.pWorkSpace[vert+i];
799 v++;
800 for ( j = 1; j < vv[0]; j++ ) {
801 for ( jj = 0; jj < nargs; jj++ ) {
802 if ( *v == arglist[jj] ) break;
803 }
804 if ( jj >= nargs ) { /* was not in the list */
805 vv[0] = vv[0]-1;
806 for ( jj = j; jj < vv[0]; jj++ ) vv[jj] = vv[jj+1];
807 }
808 v++;
809 }
810 }
811/*
812 Next we need to remove vertices that have only one object remaining.
813 Also, the object can be removed from arglist. After that we go back
814 to clean up the list.
815*/
816 for ( i = 0; i < nvert; i++ ) {
817 vv = AT.pWorkSpace[vert+i];
818 if ( vv[0] == 1 ) {
819 AT.pWorkSpace[vert+i] = AT.pWorkSpace[vert+nvert-1];
820 nvert--; i--; continue;
821 }
822 else if ( vv[0] == 2 ) {
823 for ( j = 0; j < nargs; j++ ) {
824 if ( arglist[j] == vv[1] ) break;
825 }
826 while ( j < nargs-1 ) arglist[j] = arglist[j+1];
827 nargs--;
828 AT.pWorkSpace[vert+i] = AT.pWorkSpace[vert+nvert-1];
829 nvert--; i--;
830 action = 1;
831 }
832 }
833 } while ( action );
834/*
835 Because the user can remove tadpoles rather easily as in
836 id v(?a,i?,?b,i?,?c) = v(?a,?b,?c)
837 there is no need to avoid them as potential loops.
838
839 At this point we are ready to look for loops.
840 We have
841 arglist,nargs list of eligible objects.
842 vert,nvert position in pWorkSpace for vertices and their number.
843 The vertices themselves are in the WorkSpace.
844 loop,nloop The buildup of the loop.
845*/
846 loop = AT.WorkPointer;
847 AT.WorkPointer += nargs;
848 nloop = 0;
849 numgenerated += StartLoops(BHEAD term,level,vert,nvert,arglist,nargs,loop,nloop);
850 AT.WorkPointer = oldworkpointer;
851 AT.pWorkPointer = oldpworkpointer;
852
853 if ( numgenerated == 0 && option2 != 0 ) return(Generator(BHEAD term,level));
854 return(0);
855}
856
857/*
858 #] AllLoops :
859 #[ StartLoops :
860
861 Algorithm.
862 1: have a list of vertices and objects that can be part of the loops.
863 2: starting with the first element of arglist, look for the two
864 vertices that contain this object. The first is vstart.
865 3: At the second, look at all possible objects to continue from here.
866 4: For each, find the vertex with its partner.
867 5: if the partner vertex is vstart, we have a loop. --> 9.
868 6: The partner vertex is not allowed to be a vertex that we passed
869 in the current build-up of the loop.
870 7: Put the object in loop and increase nloop.
871 8: Now sum over all allowed objects in the partner vertex. --> 4.
872
873 9: Create the function loop with the arguments from loop,nloop.
874 10: If a reverse cyclic permutation gives a smaller result, drop this
875 solution and continue with the last summing over objects at a vertex.
876 11: If not, output the result by adding the function loop
877 to the term and call Generator(term,level)
878
879 12: when all possibilities with the first element of arglist have been
880 exhausted, raise arglist and lower nargs by one. --> 2
881 13: when nargs == 1, we can stop.
882
883 The inclusion of a new object when we sum over the possibilities at a
884 vertex can best be done in a recursion, because we have no idea how
885 many nested loops we would otherwise need. Of course the recursion can
886 be emulated, but that makes the algorithm rather messy.
887
888 We have two routines: StartLoops and GenLoops.
889 StartLoops takes care of the first argument (and the summing over who
890 is the first argument), while GenLoops treats an additional vertex until
891 a loop is closed. GenLoops also sends completed loops off toe Generator.
892*/
893
894LONG StartLoops(PHEAD WORD *term,WORD level,LONG vert,WORD nvert,
895 WORD *arglist,WORD nargs,WORD *loop,WORD nloop)
896{
897 LONG numgenerated = 0;
898 WORD *v, *vv, *vstart, istart, *vpartner, ipartner, j;
899 while ( nargs > 1 ) {
900/*
901 Look for a vertex with arglist[0]. This is our starting vertex.
902 Next we look for its partner and we put arglist[0] in loop.
903*/
904 nloop = 0;
905 loop[nloop++] = arglist[0];
906 for ( istart = 0; istart < nvert; istart++ ) {
907 vstart = AT.pWorkSpace[vert+istart];
908 v = vstart+1; vv = vstart + *vstart;
909 do {
910 if ( *v == arglist[0] ) goto havestart;
911 v++;
912 } while ( v < vv );
913 }
914/*
915 If we come here, we have a problem.
916*/
917/* INTERNAL_ERROR_EXCL_START */
918 MesPrint("!>Internal error in StartLoops. Object not found.");
919 Terminate(-1);
920 return(-1);
921/* INTERNAL_ERROR_EXCL_STOP */
922havestart:
923 AT.pWorkSpace[vert+nvert] = vstart;
924/*
925 Check for tadpole.
926*/
927 v++;
928 while ( v < vv ) {
929 if ( *v == arglist[0] ) { /* tadpole */
930 LoopOutput(BHEAD term,level,loop,nloop);
931 numgenerated++;
932 goto nextarg;
933 }
934 v++;
935 }
936/*
937 Now the partner vertex.
938*/
939 for ( ipartner = istart+1; ipartner < nvert; ipartner++ ) {
940 vpartner = AT.pWorkSpace[vert+ipartner];
941 vv = vpartner+*vpartner; v = vpartner+1;
942 do {
943 if ( *v == arglist[0] ) goto havepartner;
944 v++;
945 } while ( v < vv );
946 }
947 return(numgenerated);
948havepartner:
949 AT.pWorkSpace[vert+nvert+1] = vpartner;
950/*
951 Now we run through all other possibilities at vpartner.
952 vv is still OK.
953*/
954 v = vpartner+1;
955 while ( v < vv ) {
956 if ( *v != arglist[0] ) {
957 for ( j = 1; j < nargs; j++ ) {
958 if ( *v == arglist[j] ) {
959 loop[nloop++] = *v;
960 numgenerated += GenLoops(BHEAD term,level,vert,nvert,
961 arglist,nargs,loop,nloop);
962 nloop--;
963 break;
964 }
965 }
966 }
967 v++;
968 }
969nextarg:
970 arglist++; nargs--;
971 }
972 return(numgenerated);
973}
974
975/*
976 #] StartLoops :
977 #[ GenLoops :
978
979 We enter with an open line in loop[nloop-1]
980*/
981
982LONG GenLoops(PHEAD WORD *term,WORD level,LONG vert,WORD nvert,
983 WORD *arglist,WORD nargs,WORD *loop,WORD nloop)
984{
985 LONG numgenerated = 0;
986 WORD *vstart, *v, *vv, i, j, *vpartner;
987/*
988 Start with checking whether the partner is in vstart (=vert[nvert])
989*/
990 vstart = AT.pWorkSpace[nvert];
991 vv = vstart + *vstart; v = vstart+1;
992 while ( v < vv ) {
993 if ( *v == loop[nloop-1] ) {
994/*
995 This closes the loop.
996 Now we can output it.
997*/
998 LoopOutput(BHEAD term,level,loop,nloop);
999 numgenerated++;
1000 return(numgenerated);
1001 }
1002 v++;
1003 }
1004/*
1005 Start with finding the partner.
1006*/
1007 for ( i = 0; i < nvert; i++ ) {
1008 vpartner = AT.pWorkSpace[vert+i];
1009 if ( vpartner == vstart ) continue;
1010 for ( j = 0; j < nloop; j++ ) {
1011 if ( vpartner == AT.pWorkSpace[vert+nvert+j] ) break;
1012 }
1013 if ( j < nloop ) continue;
1014 v = vpartner+1; vv = vpartner + *vpartner;
1015 while ( v < vv ) {
1016 if ( *v == loop[nloop-1] ) {
1017/*
1018 Found the partner.
1019 Now we can sum over the remaining permitted arguments of this vertex.
1020*/
1021 v = vpartner + 1;
1022 while ( v < vv ) {
1023/*
1024 *v should be in arglist.
1025*/
1026 for ( j = 0; j < nargs; j++ ) {
1027 if ( *v == arglist[j] ) break;
1028 }
1029 if ( j >= nargs ) { v++; continue; }
1030/*
1031 *v should not be in loops.
1032*/
1033 for ( j = 0; j < nloop; j++ ) {
1034 if ( *v == loop[j] ) break;
1035 }
1036 if ( j >= nloop ) {
1037 AT.pWorkSpace[vert+nvert+nloop] = vpartner;
1038 loop[nloop++] = *v;
1039 numgenerated += GenLoops(BHEAD term,level,vert,nvert,
1040 arglist,nargs,loop,nloop);
1041 nloop--;
1042 }
1043 v++;
1044 }
1045 return(numgenerated);
1046 }
1047 v++;
1048 }
1049 }
1050/*
1051 The partner is in a vertex we have finished before.
1052*/
1053 return(numgenerated);
1054}
1055
1056/*
1057 #] GenLoops :
1058 #[ LoopOutput :
1059*/
1060
1061void LoopOutput(PHEAD WORD *term, WORD level, WORD *loop, WORD nloop)
1062{
1063 CBUF *C = cbuf+AM.rbufnum;
1064 WORD loopfun = C->lhs[level][3]; /* The output function */
1065 WORD option1 = C->lhs[level][4]; /* type of argument */
1066 WORD *tstop, *tstop1, *t, *tt;
1067 WORD *outterm, *loop1;
1068 WORD i;
1069 tstop1 = term+*term; tstop = tstop1 - ABS(tstop1[-1]);
1070/*
1071 Construct the rcycle symmetrized version of loop.
1072*/
1073 if ( nloop > 2 ) {
1074 loop1 = AT.WorkPointer;
1075 loop1[0] = loop[0];
1076 for ( i = 1; i < nloop; i++ ) { loop1[i] = loop[nloop-i]; }
1077 if ( loop1[1] < loop[1] ) {
1078 AT.WorkPointer += nloop;
1079 loop = loop1;
1080 }
1081 }
1082 outterm = AT.WorkPointer;
1083 tt = outterm; t = term;
1084 while ( t < tstop ) *tt++ = *t++;
1085 *tt++ = loopfun;
1086 if ( functions[loopfun-FUNCTION].spec == TENSORFUNCTION ) {
1087 *tt++ = FUNHEAD+nloop;
1088 FILLFUN(tt)
1089 if ( option1 == -VECTOR ) {
1090 for ( i = 0; i < nloop; i++ ) *tt++ = loop[i]+AM.OffsetVector;
1091 }
1092 else {
1093 for ( i = 0; i < nloop; i++ ) *tt++ = loop[i];
1094 }
1095 }
1096 else {
1097 *tt++ = FUNHEAD+nloop*2;
1098 FILLFUN(tt)
1099 for ( i = 0; i< nloop; i++ ) {
1100 *tt++ = option1;
1101 if ( option1 == -VECTOR ) *tt++ = loop[i] + AM.OffsetVector;
1102 else *tt++ = loop[i];
1103 }
1104 }
1105 while ( t < tstop1 ) *tt++ = *t++;
1106 *outterm = tt - outterm;
1107 AT.WorkPointer = tt;
1108 if ( Generator(BHEAD outterm,level) ) {
1109 MesCall("LoopOutput");
1110 Terminate(-1);
1111 }
1112 AT.WorkPointer = outterm;
1113}
1114
1115/*
1116 #] LoopOutput :
1117 #[ AllPaths :
1118
1119 This routine has many similarities with AllLoops.
1120 In AllLoops we have startingpoint and endpoint at the same vertex.
1121 In AllPaths the startingpoint and the endpoints are different and
1122 given in advance. This endfun can occur only twice.
1123*/
1124
1125int AllPaths(PHEAD WORD *term,WORD level)
1126{
1127 CBUF *C = cbuf+AM.rbufnum;
1128 WORD endcode = C->lhs[level][2]; /* The endpoint function */
1129 WORD vcode = C->lhs[level][3]; /* The intermediate function */
1130 WORD option1 = C->lhs[level][5]; /* type of argument */
1131 WORD option2 = C->lhs[level][6]; /* what to do when no loop */
1132 WORD *t, *tstop, *tend1, *tend2, *tstart, *tend, numend, nvert, npass;
1133 WORD *tfrom, *to, *tos, *from;
1134 WORD *arglist, nargs, *path, npath, *a, *a1, *a2, *a3;
1135 WORD i, j, jj, *v, *vv, action;
1136 LONG vert,vert1; /* ,vert2; */
1137 WORD numgenerated = 0;
1138 WORD *oldworkpointer = AT.WorkPointer;
1139 LONG oldpworkpointer = AT.pWorkPointer;
1140/*
1141 Search for the two occurrences of endcode.
1142*/
1143 tstop = term+*term; tstop -= ABS(tstop[-1]);
1144 t = term + 1;
1145 while ( t < tstop && *t != endcode ) t += t[1];
1146 if ( t == tstop ) { /* no endpoints */
1147 if ( option2 == 0 ) return(0);
1148 else return(Generator(BHEAD term,level));
1149 }
1150 tend1 = t;
1151 numend = 0;
1152 while ( t < tstop && *t == endcode ) { tend2 = t; t += t[1]; numend++; }
1153 if ( numend != 2 ) { /* not 2 endpoints -> no path */
1154 if ( option2 == 0 ) return(0);
1155 else return(Generator(BHEAD term,level));
1156 }
1157/*
1158 Search for the intermediate functions.
1159*/
1160 t = term + 1;
1161 nvert = 0;
1162 while ( t < tstop && *t != vcode ) t += t[1];
1163 tstart = t;
1164 while ( t < tstop && *t == vcode ) {
1165 t += t[1]; nvert++;
1166 }
1167 tend = t;
1168/*
1169 Next we copy the relevant content of the various functions into the
1170 WorkSpace. We keep pointers to the functions in pWorkSpace.
1171 First the two endpoints and then the intermediate ones.
1172*/
1173 WantAddPointers((2*nvert+8));
1174 vert = AT.pWorkPointer+2;
1175 vert1 = AT.pWorkPointer;
1176/* vert2 = AT.pWorkPointer+1; */
1177 AT.pWorkPointer += 2*nvert+8;
1178 nvert = 0;
1179/*
1180 Copy the endpoints.
1181*/
1182 to = AT.WorkPointer;
1183
1184 if ( functions[endcode-FUNCTION].spec == TENSORFUNCTION ) {
1185 from = tend1; npass = 0;
1186redo1:
1187 tos = to++;
1188 tfrom = from+from[1];
1189 from += FUNHEAD;
1190 while ( from < tfrom ) {
1191 if ( option1 == -INDEX && *from >= 0 ) {
1192 *to++ = *from++;
1193 }
1194 else if ( option1 == -VECTOR && *from < 0 ) {
1195 *to++ = *from++ - AM.OffsetVector;
1196 }
1197 else {
1198 from++;
1199 }
1200 }
1201 *tos = to-tos;
1202 if ( *tos < 2 ) to = tos;
1203 else AT.pWorkSpace[vert1+npass] = tos;
1204 npass++;
1205 if ( from == tend2 ) goto redo1;
1206 }
1207 else if ( functions[endcode-FUNCTION].spec == 0 ) {
1208 from = tend1;
1209 npass = 0;
1210redo2:
1211 tos = to++;
1212 tfrom = from + from[1];
1213 from += FUNHEAD;
1214 while ( from < tfrom ) {
1215 if ( option1 == -VECTOR
1216 && ( from[0] == -INDEX || from[0] == -VECTOR )
1217 && from[1] < 0 ) {
1218 from++;
1219 *to++ = *from++ - AM.OffsetVector;
1220 }
1221 else if ( option1 == *from ) {
1222 from++; *to++ = *from++;
1223 }
1224 else {
1225 NEXTARG(from);
1226 }
1227 }
1228 *tos = to-tos;
1229 if ( *tos < 2 ) to = tos;
1230 else AT.pWorkSpace[vert1+npass] = tos;
1231 npass++;
1232 if ( from == tend2 ) goto redo2;
1233 }
1234 else {
1235 AT.WorkPointer = oldworkpointer;
1236 AT.pWorkPointer = oldpworkpointer;
1237 if ( option2 == 0 ) return(0);
1238 return(Generator(BHEAD term,level));
1239 }
1240/*
1241 And now the intermediate functions
1242*/
1243 from = tstart;
1244 if ( functions[vcode-FUNCTION].spec == TENSORFUNCTION ) {
1245 from = tstart;
1246 while ( from < tend ) {
1247 tos = to++;
1248 tfrom = from+from[1];
1249 from += FUNHEAD;
1250 while ( from < tfrom ) {
1251 if ( option1 == -INDEX && *from >= 0 ) {
1252 *to++ = *from++;
1253 }
1254 else if ( option1 == -VECTOR && *from < 0 ) {
1255 *to++ = *from++ - AM.OffsetVector;
1256 }
1257 else {
1258 from++;
1259 }
1260 }
1261 *tos = to-tos;
1262 if ( *tos < 3 ) to = tos;
1263 else AT.pWorkSpace[vert+nvert++] = tos;
1264 }
1265 }
1266 else if ( functions[vcode-FUNCTION].spec == 0 ) {
1267 from = tstart;
1268 while ( from < tend ) {
1269 tos = to++;
1270 tfrom = from + from[1];
1271 from += FUNHEAD;
1272 while ( from < tfrom ) {
1273 if ( option1 == -VECTOR
1274 && ( from[0] == -INDEX || from[0] == -VECTOR )
1275 && from[1] < 0 ) {
1276 from++;
1277 *to++ = *from++ - AM.OffsetVector;
1278 }
1279 else if ( option1 == *from ) {
1280 from++; *to++ = *from++;
1281 }
1282 else {
1283 NEXTARG(from);
1284 }
1285 }
1286 *tos = to-tos;
1287 if ( *tos < 3 ) to = tos;
1288 else AT.pWorkSpace[vert+nvert++] = tos;
1289 }
1290 }
1291 else {
1292 AT.WorkPointer = oldworkpointer;
1293 AT.pWorkPointer = oldpworkpointer;
1294 if ( option2 == 0 ) return(0);
1295 return(Generator(BHEAD term,level));
1296 }
1297 AT.WorkPointer = to;
1298/*
1299 Now make a list of all relevant arguments.
1300*/
1301 a = arglist = AT.WorkPointer;
1302 for ( i = -2; i < nvert; i++ ) {
1303 v = AT.pWorkSpace[vert+i];
1304 j = *v-1; v++;
1305 NCOPY(a,v,j);
1306 }
1307 nargs = a-arglist;
1308 AT.WorkPointer = a;
1309/*
1310 Now sort the list. Bubble type sort.
1311*/
1312 a1 = arglist; a2 = a1+1;
1313 while ( a2 < a ) {
1314 a3 = a2;
1315 while ( a3 > a1 && a3[-1] > a3[0] ) {
1316 EXCH(*a3,a3[-1]);
1317 a3--;
1318 }
1319 a2++;
1320 }
1321/*
1322 Now remove the elements from the list that do not occur exactly twice.
1323*/
1324 a1 = arglist; a2 = a1; a3 = a1+nargs;
1325 while ( a2 < a3 ) {
1326 if ( a2+1 == a3 ) { break; }
1327 else if ( a2+2 == a3 ) {
1328 if ( a2[0] == a2[1] ) { *a1++ = a2[0]; }
1329 break;
1330 }
1331 else {
1332 if ( a2[0] != a2[1] ) { a2++; }
1333 else if ( a2[0] != a2[2] ) {
1334 *a1++ = a2[0]; a2 += 2;
1335 }
1336 else {
1337 a2++; while ( a2 < a3 && a2[-1] == a2[0] ) a2++;
1338 }
1339 }
1340 }
1341 nargs = a1-arglist;
1342/*
1343 Now we need to redo the list of vertices and remove the elements
1344 that are not in our arglist.
1345*/
1346 do {
1347 action = 0;
1348 for ( i = -2; i < nvert; i++ ) {
1349 vv = v = AT.pWorkSpace[vert+i];
1350 v++;
1351 for ( j = 1; j < vv[0]; j++ ) {
1352 for ( jj = 0; jj < nargs; jj++ ) {
1353 if ( *v == arglist[jj] ) break;
1354 }
1355 if ( jj >= nargs ) { /* was not in the list */
1356 vv[0] = vv[0]-1;
1357 for ( jj = j; jj < vv[0]; jj++ ) vv[jj] = vv[jj+1];
1358 }
1359 v++;
1360 }
1361 }
1362/*
1363 Next we need to remove vertices that have only one object remaining.
1364 Also, the object can be removed from arglist. After that we go back
1365 to clean up the list.
1366*/
1367 for ( i = -2; i < nvert; i++ ) {
1368 vv = AT.pWorkSpace[vert+i];
1369 if ( vv[0] == 1 ) {
1370 AT.pWorkSpace[vert+i] = AT.pWorkSpace[vert+nvert-1];
1371 nvert--; i--; continue;
1372 }
1373 else if ( vv[0] == 2 && i >= 0 ) {
1374 for ( j = 0; j < nargs; j++ ) {
1375 if ( arglist[j] == vv[1] ) break;
1376 }
1377 while ( j < nargs-1 ) arglist[j] = arglist[j+1];
1378 nargs--;
1379 AT.pWorkSpace[vert+i] = AT.pWorkSpace[vert+nvert-1];
1380 nvert--; i--;
1381 action = 1;
1382 }
1383 }
1384 } while ( action );
1385/*
1386 Now we have a clean list of connections and we can start building
1387 the path. We start with summing over all objects in vert1.
1388*/
1389 path = AT.WorkPointer; npath = 0;
1390 AT.WorkPointer += nvert+8;
1391
1392 t = AT.pWorkSpace[vert1];
1393 if ( *t >= 2 ) {
1394 for ( i = 1; i < *t; i++ ) {
1395 AT.pWorkSpace[vert+nvert] = t;
1396 path[npath++] = t[i];
1397 numgenerated += GenPaths(BHEAD term,level,vert,nvert,arglist,nargs,path,npath);
1398 npath--;
1399 }
1400 }
1401 AT.WorkPointer = oldworkpointer;
1402 AT.pWorkPointer = oldpworkpointer;
1403 if ( numgenerated == 0 && option2 != 0 ) return(Generator(BHEAD term,level));
1404 return(0);
1405}
1406
1407/*
1408 #] AllPaths :
1409 #[ GenPaths :
1410
1411 We enter with an open line in path[npath-1]
1412 We traverse though the vertices until we reach vert-1, the endpoint.
1413*/
1414
1415LONG GenPaths(PHEAD WORD *term, WORD level, LONG vert, WORD nvert,
1416 WORD *arglist, WORD nargs, WORD *path, WORD npath)
1417{
1418 LONG numgenerated = 0;
1419 WORD *t, *vpartner, *v, *vv;
1420 WORD i, j;
1421/*
1422 Check whether path[npath-1] is part of the endpoint.
1423*/
1424 t = AT.pWorkSpace[vert-1];
1425 for ( i = 1; i < *t; i++ ) {
1426 if ( t[i] == path[npath-1] ) { /* Got a path! */
1427 PathOutput(BHEAD term,level,path,npath);
1428 numgenerated++;
1429 return(numgenerated);
1430 }
1431 }
1432/*
1433 Start with finding the partner.
1434*/
1435 for ( i = 0; i < nvert; i++ ) {
1436 vpartner = AT.pWorkSpace[vert+i];
1437 for ( j = 0; j < npath; j++ ) {
1438 if ( vpartner == AT.pWorkSpace[vert+nvert+j] ) break;
1439 }
1440 if ( j < npath ) continue;
1441 v = vpartner+1; vv = vpartner + *vpartner;
1442 while ( v < vv ) {
1443 if ( *v == path[npath-1] ) {
1444/*
1445 Found the partner.
1446 Now we can sum over the remaining permitted arguments of this vertex.
1447*/
1448 v = vpartner + 1;
1449 while ( v < vv ) {
1450/*
1451 *v should be in arglist.
1452*/
1453 for ( j = 0; j < nargs; j++ ) {
1454 if ( *v == arglist[j] ) break;
1455 }
1456 if ( j >= nargs ) { v++; continue; }
1457/*
1458 *v should not be in path.
1459*/
1460 for ( j = 0; j < npath; j++ ) {
1461 if ( *v == path[j] ) break;
1462 }
1463 if ( j >= npath ) {
1464 AT.pWorkSpace[vert+nvert+npath] = vpartner;
1465 path[npath++] = *v;
1466 numgenerated += GenPaths(BHEAD term,level,vert,nvert,
1467 arglist,nargs,path,npath);
1468 npath--;
1469 }
1470 v++;
1471 }
1472 return(numgenerated);
1473 }
1474 v++;
1475 }
1476 }
1477/*
1478 The partner is in a vertex we have finished before.
1479*/
1480 return(numgenerated);
1481}
1482
1483/*
1484 #] GenPaths :
1485 #[ PathOutput :
1486*/
1487
1488void PathOutput(PHEAD WORD *term, WORD level, WORD *path, WORD npath)
1489{
1490 CBUF *C = cbuf+AM.rbufnum;
1491 WORD pathfun = C->lhs[level][4]; /* The output function */
1492 WORD option1 = C->lhs[level][5]; /* type of argument */
1493 WORD *tstop, *tstop1, *t, *tt;
1494 WORD *outterm;
1495 WORD i;
1496 tstop1 = term+*term; tstop = tstop1 - ABS(tstop1[-1]);
1497 outterm = AT.WorkPointer;
1498 tt = outterm; t = term;
1499 while ( t < tstop ) *tt++ = *t++;
1500 *tt++ = pathfun;
1501 if ( functions[pathfun-FUNCTION].spec == TENSORFUNCTION ) {
1502 *tt++ = FUNHEAD+npath;
1503 FILLFUN(tt)
1504 if ( option1 == -VECTOR ) {
1505 for ( i = 0; i < npath; i++ ) *tt++ = path[i]+AM.OffsetVector;
1506 }
1507 else {
1508 for ( i = 0; i < npath; i++ ) *tt++ = path[i];
1509 }
1510 }
1511 else {
1512 *tt++ = FUNHEAD+npath*2;
1513 FILLFUN(tt)
1514 for ( i = 0; i< npath; i++ ) {
1515 *tt++ = option1;
1516 if ( option1 == -VECTOR ) *tt++ = path[i] + AM.OffsetVector;
1517 else *tt++ = path[i];
1518 }
1519 }
1520 while ( t < tstop1 ) *tt++ = *t++;
1521 *outterm = tt - outterm;
1522 AT.WorkPointer = tt;
1523 if ( Generator(BHEAD outterm,level) ) {
1524 MesCall("PathOutput");
1525 Terminate(-1);
1526 }
1527 AT.WorkPointer = outterm;
1528}
1529
1530
1531
1532/*
1533 #] PathOutput :
1534 #[ AllOnePI :
1535
1536 We assume a graph that has loops and is OnePI.
1537 This routine initializes the recursion that creates all onePI subgraphs.
1538
1539 Algorithm:
1540 a: have a routine that removes all bridges: RemoveBridges
1541 b: output an empty diagram.
1542 c: output the diagram itself
1543 d: cut a line, followed by removing all bridges.
1544 e: if the diagram still has lines, go to c, else go to the next line in d.
1545 In order to avoid factorial blowup, we need to prune the tree as fast as possible.
1546
1547 1: list of lines to be cut.
1548 2: once we have tried a line and removed its bridges, we do not have to
1549 try this line deeper in the tree, neither its bridges.
1550
1551
1552 1 2 3
1553 cut 1. x
1554 cut 2 -> bridge 3 x
1555 cut 2. x
1556 cut 3 -> bridge 1 ???? mag niet
1557 cut 3. x
1558 Hence 1,2,3,4,5,6,7
1559 1 still to go 2,3,4,5,6,7
1560 2 still to go 3,4,5,6,7
1561 3 still to go 4,5,6,7
1562 etc.
1563 bridges: if x is bridge, we take x from the list_to_go.
1564 If we have a bridge that is a number lower than the list_to_go we skip this possibility.
1565 We reach the end when the list_to_go is empty.
1566*/
1567
1568WORD AllOnePI(WORD *term,WORD level)
1569{
1570 CBUF *C = cbuf+AM.rbufnum;
1571 WORD vcode = C->lhs[level][2]; /* The vertex function */
1572 WORD option1 = C->lhs[level][4]; /* type of argument */
1573/*
1574 First we have to collect all relevant information about the diagram.
1575 We should start with removing bridges to make the real starting point onePI.
1576*/
1577 DUMMYUSE(term)
1578 DUMMYUSE(vcode)
1579 DUMMYUSE(option1)
1580 return(0);
1581}
1582
1583/*
1584 #] AllOnePI :
1585 #[ RemoveBridges :
1586*/
1587
1588int RemoveBridges(void)
1589{
1590 return(0);
1591}
1592
1593/*
1594 #] RemoveBridges :
1595 #[ TakeOneLine :
1596*/
1597
1598int TakeOneLine(WORD*term,WORD level)
1599{
1600 DUMMYUSE(term)
1601 DUMMYUSE(level)
1602 return(0);
1603}
1604
1605/*
1606 #] TakeOneLine :
1607 #[ OutputOnePI :
1608*/
1609
1610int OutputOnePI(PHEAD WORD *term,WORD level)
1611{
1612 return(Generator(BHEAD term,level));
1613}
1614
1615/*
1616 #] OutputOnePI :
1617*/
1618
Generator
int Generator(PHEAD WORD *, WORD)
Definition proces.c:3275
CbUf::lhs
WORD ** lhs
Definition structs.h:974
Generated by doxygen 1.9.8