FORM v5.0.0-35-g6318119
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 MesPrint("Internal error in StartLoops. Object not found.");
918 Terminate(-1);
919 return(-1);
920havestart:
921 AT.pWorkSpace[vert+nvert] = vstart;
922/*
923 Check for tadpole.
924*/
925 v++;
926 while ( v < vv ) {
927 if ( *v == arglist[0] ) { /* tadpole */
928 LoopOutput(BHEAD term,level,loop,nloop);
929 numgenerated++;
930 goto nextarg;
931 }
932 v++;
933 }
934/*
935 Now the partner vertex.
936*/
937 for ( ipartner = istart+1; ipartner < nvert; ipartner++ ) {
938 vpartner = AT.pWorkSpace[vert+ipartner];
939 vv = vpartner+*vpartner; v = vpartner+1;
940 do {
941 if ( *v == arglist[0] ) goto havepartner;
942 v++;
943 } while ( v < vv );
944 }
945 return(numgenerated);
946havepartner:
947 AT.pWorkSpace[vert+nvert+1] = vpartner;
948/*
949 Now we run through all other possibilities at vpartner.
950 vv is still OK.
951*/
952 v = vpartner+1;
953 while ( v < vv ) {
954 if ( *v != arglist[0] ) {
955 for ( j = 1; j < nargs; j++ ) {
956 if ( *v == arglist[j] ) {
957 loop[nloop++] = *v;
958 numgenerated += GenLoops(BHEAD term,level,vert,nvert,
959 arglist,nargs,loop,nloop);
960 nloop--;
961 break;
962 }
963 }
964 }
965 v++;
966 }
967nextarg:
968 arglist++; nargs--;
969 }
970 return(numgenerated);
971}
972
973/*
974 #] StartLoops :
975 #[ GenLoops :
976
977 We enter with an open line in loop[nloop-1]
978*/
979
980LONG GenLoops(PHEAD WORD *term,WORD level,LONG vert,WORD nvert,
981 WORD *arglist,WORD nargs,WORD *loop,WORD nloop)
982{
983 LONG numgenerated = 0;
984 WORD *vstart, *v, *vv, i, j, *vpartner;
985/*
986 Start with checking whether the partner is in vstart (=vert[nvert])
987*/
988 vstart = AT.pWorkSpace[nvert];
989 vv = vstart + *vstart; v = vstart+1;
990 while ( v < vv ) {
991 if ( *v == loop[nloop-1] ) {
992/*
993 This closes the loop.
994 Now we can output it.
995*/
996 LoopOutput(BHEAD term,level,loop,nloop);
997 numgenerated++;
998 return(numgenerated);
999 }
1000 v++;
1001 }
1002/*
1003 Start with finding the partner.
1004*/
1005 for ( i = 0; i < nvert; i++ ) {
1006 vpartner = AT.pWorkSpace[vert+i];
1007 if ( vpartner == vstart ) continue;
1008 for ( j = 0; j < nloop; j++ ) {
1009 if ( vpartner == AT.pWorkSpace[vert+nvert+j] ) break;
1010 }
1011 if ( j < nloop ) continue;
1012 v = vpartner+1; vv = vpartner + *vpartner;
1013 while ( v < vv ) {
1014 if ( *v == loop[nloop-1] ) {
1015/*
1016 Found the partner.
1017 Now we can sum over the remaining permitted arguments of this vertex.
1018*/
1019 v = vpartner + 1;
1020 while ( v < vv ) {
1021/*
1022 *v should be in arglist.
1023*/
1024 for ( j = 0; j < nargs; j++ ) {
1025 if ( *v == arglist[j] ) break;
1026 }
1027 if ( j >= nargs ) { v++; continue; }
1028/*
1029 *v should not be in loops.
1030*/
1031 for ( j = 0; j < nloop; j++ ) {
1032 if ( *v == loop[j] ) break;
1033 }
1034 if ( j >= nloop ) {
1035 AT.pWorkSpace[vert+nvert+nloop] = vpartner;
1036 loop[nloop++] = *v;
1037 numgenerated += GenLoops(BHEAD term,level,vert,nvert,
1038 arglist,nargs,loop,nloop);
1039 nloop--;
1040 }
1041 v++;
1042 }
1043 return(numgenerated);
1044 }
1045 v++;
1046 }
1047 }
1048/*
1049 The partner is in a vertex we have finished before.
1050*/
1051 return(numgenerated);
1052}
1053
1054/*
1055 #] GenLoops :
1056 #[ LoopOutput :
1057*/
1058
1059void LoopOutput(PHEAD WORD *term, WORD level, WORD *loop, WORD nloop)
1060{
1061 CBUF *C = cbuf+AM.rbufnum;
1062 WORD loopfun = C->lhs[level][3]; /* The output function */
1063 WORD option1 = C->lhs[level][4]; /* type of argument */
1064 WORD *tstop, *tstop1, *t, *tt;
1065 WORD *outterm, *loop1;
1066 WORD i;
1067 tstop1 = term+*term; tstop = tstop1 - ABS(tstop1[-1]);
1068/*
1069 Construct the rcycle symmetrized version of loop.
1070*/
1071 if ( nloop > 2 ) {
1072 loop1 = AT.WorkPointer;
1073 loop1[0] = loop[0];
1074 for ( i = 1; i < nloop; i++ ) { loop1[i] = loop[nloop-i]; }
1075 if ( loop1[1] < loop[1] ) {
1076 AT.WorkPointer += nloop;
1077 loop = loop1;
1078 }
1079 }
1080 outterm = AT.WorkPointer;
1081 tt = outterm; t = term;
1082 while ( t < tstop ) *tt++ = *t++;
1083 *tt++ = loopfun;
1084 if ( functions[loopfun-FUNCTION].spec == TENSORFUNCTION ) {
1085 *tt++ = FUNHEAD+nloop;
1086 FILLFUN(tt)
1087 if ( option1 == -VECTOR ) {
1088 for ( i = 0; i < nloop; i++ ) *tt++ = loop[i]+AM.OffsetVector;
1089 }
1090 else {
1091 for ( i = 0; i < nloop; i++ ) *tt++ = loop[i];
1092 }
1093 }
1094 else {
1095 *tt++ = FUNHEAD+nloop*2;
1096 FILLFUN(tt)
1097 for ( i = 0; i< nloop; i++ ) {
1098 *tt++ = option1;
1099 if ( option1 == -VECTOR ) *tt++ = loop[i] + AM.OffsetVector;
1100 else *tt++ = loop[i];
1101 }
1102 }
1103 while ( t < tstop1 ) *tt++ = *t++;
1104 *outterm = tt - outterm;
1105 AT.WorkPointer = tt;
1106 if ( Generator(BHEAD outterm,level) ) {
1107 MesCall("LoopOutput");
1108 Terminate(-1);
1109 }
1110 AT.WorkPointer = outterm;
1111}
1112
1113/*
1114 #] LoopOutput :
1115 #[ AllPaths :
1116
1117 This routine has many similarities with AllLoops.
1118 In AllLoops we have startingpoint and endpoint at the same vertex.
1119 In AllPaths the startingpoint and the endpoints are different and
1120 given in advance. This endfun can occur only twice.
1121*/
1122
1123int AllPaths(PHEAD WORD *term,WORD level)
1124{
1125 CBUF *C = cbuf+AM.rbufnum;
1126 WORD endcode = C->lhs[level][2]; /* The endpoint function */
1127 WORD vcode = C->lhs[level][3]; /* The intermediate function */
1128 WORD option1 = C->lhs[level][5]; /* type of argument */
1129 WORD option2 = C->lhs[level][6]; /* what to do when no loop */
1130 WORD *t, *tstop, *tend1, *tend2, *tstart, *tend, numend, nvert, npass;
1131 WORD *tfrom, *to, *tos, *from;
1132 WORD *arglist, nargs, *path, npath, *a, *a1, *a2, *a3;
1133 WORD i, j, jj, *v, *vv, action;
1134 LONG vert,vert1; /* ,vert2; */
1135 WORD numgenerated = 0;
1136 WORD *oldworkpointer = AT.WorkPointer;
1137 LONG oldpworkpointer = AT.pWorkPointer;
1138/*
1139 Search for the two occurrences of endcode.
1140*/
1141 tstop = term+*term; tstop -= ABS(tstop[-1]);
1142 t = term + 1;
1143 while ( t < tstop && *t != endcode ) t += t[1];
1144 if ( t == tstop ) { /* no endpoints */
1145 if ( option2 == 0 ) return(0);
1146 else return(Generator(BHEAD term,level));
1147 }
1148 tend1 = t;
1149 numend = 0;
1150 while ( t < tstop && *t == endcode ) { tend2 = t; t += t[1]; numend++; }
1151 if ( numend != 2 ) { /* not 2 endpoints -> no path */
1152 if ( option2 == 0 ) return(0);
1153 else return(Generator(BHEAD term,level));
1154 }
1155/*
1156 Search for the intermediate functions.
1157*/
1158 t = term + 1;
1159 nvert = 0;
1160 while ( t < tstop && *t != vcode ) t += t[1];
1161 tstart = t;
1162 while ( t < tstop && *t == vcode ) {
1163 t += t[1]; nvert++;
1164 }
1165 tend = t;
1166/*
1167 Next we copy the relevant content of the various functions into the
1168 WorkSpace. We keep pointers to the functions in pWorkSpace.
1169 First the two endpoints and then the intermediate ones.
1170*/
1171 WantAddPointers((2*nvert+8));
1172 vert = AT.pWorkPointer+2;
1173 vert1 = AT.pWorkPointer;
1174/* vert2 = AT.pWorkPointer+1; */
1175 AT.pWorkPointer += 2*nvert+8;
1176 nvert = 0;
1177/*
1178 Copy the endpoints.
1179*/
1180 to = AT.WorkPointer;
1181
1182 if ( functions[endcode-FUNCTION].spec == TENSORFUNCTION ) {
1183 from = tend1; npass = 0;
1184redo1:
1185 tos = to++;
1186 tfrom = from+from[1];
1187 from += FUNHEAD;
1188 while ( from < tfrom ) {
1189 if ( option1 == -INDEX && *from >= 0 ) {
1190 *to++ = *from++;
1191 }
1192 else if ( option1 == -VECTOR && *from < 0 ) {
1193 *to++ = *from++ - AM.OffsetVector;
1194 }
1195 else {
1196 from++;
1197 }
1198 }
1199 *tos = to-tos;
1200 if ( *tos < 2 ) to = tos;
1201 else AT.pWorkSpace[vert1+npass] = tos;
1202 npass++;
1203 if ( from == tend2 ) goto redo1;
1204 }
1205 else if ( functions[endcode-FUNCTION].spec == 0 ) {
1206 from = tend1;
1207 npass = 0;
1208redo2:
1209 tos = to++;
1210 tfrom = from + from[1];
1211 from += FUNHEAD;
1212 while ( from < tfrom ) {
1213 if ( option1 == -VECTOR
1214 && ( from[0] == -INDEX || from[0] == -VECTOR )
1215 && from[1] < 0 ) {
1216 from++;
1217 *to++ = *from++ - AM.OffsetVector;
1218 }
1219 else if ( option1 == *from ) {
1220 from++; *to++ = *from++;
1221 }
1222 else {
1223 NEXTARG(from);
1224 }
1225 }
1226 *tos = to-tos;
1227 if ( *tos < 2 ) to = tos;
1228 else AT.pWorkSpace[vert1+npass] = tos;
1229 npass++;
1230 if ( from == tend2 ) goto redo2;
1231 }
1232 else {
1233 AT.WorkPointer = oldworkpointer;
1234 AT.pWorkPointer = oldpworkpointer;
1235 if ( option2 == 0 ) return(0);
1236 return(Generator(BHEAD term,level));
1237 }
1238/*
1239 And now the intermediate functions
1240*/
1241 from = tstart;
1242 if ( functions[vcode-FUNCTION].spec == TENSORFUNCTION ) {
1243 from = tstart;
1244 while ( from < tend ) {
1245 tos = to++;
1246 tfrom = from+from[1];
1247 from += FUNHEAD;
1248 while ( from < tfrom ) {
1249 if ( option1 == -INDEX && *from >= 0 ) {
1250 *to++ = *from++;
1251 }
1252 else if ( option1 == -VECTOR && *from < 0 ) {
1253 *to++ = *from++ - AM.OffsetVector;
1254 }
1255 else {
1256 from++;
1257 }
1258 }
1259 *tos = to-tos;
1260 if ( *tos < 3 ) to = tos;
1261 else AT.pWorkSpace[vert+nvert++] = tos;
1262 }
1263 }
1264 else if ( functions[vcode-FUNCTION].spec == 0 ) {
1265 from = tstart;
1266 while ( from < tend ) {
1267 tos = to++;
1268 tfrom = from + from[1];
1269 from += FUNHEAD;
1270 while ( from < tfrom ) {
1271 if ( option1 == -VECTOR
1272 && ( from[0] == -INDEX || from[0] == -VECTOR )
1273 && from[1] < 0 ) {
1274 from++;
1275 *to++ = *from++ - AM.OffsetVector;
1276 }
1277 else if ( option1 == *from ) {
1278 from++; *to++ = *from++;
1279 }
1280 else {
1281 NEXTARG(from);
1282 }
1283 }
1284 *tos = to-tos;
1285 if ( *tos < 3 ) to = tos;
1286 else AT.pWorkSpace[vert+nvert++] = tos;
1287 }
1288 }
1289 else {
1290 AT.WorkPointer = oldworkpointer;
1291 AT.pWorkPointer = oldpworkpointer;
1292 if ( option2 == 0 ) return(0);
1293 return(Generator(BHEAD term,level));
1294 }
1295 AT.WorkPointer = to;
1296/*
1297 Now make a list of all relevant arguments.
1298*/
1299 a = arglist = AT.WorkPointer;
1300 for ( i = -2; i < nvert; i++ ) {
1301 v = AT.pWorkSpace[vert+i];
1302 j = *v-1; v++;
1303 NCOPY(a,v,j);
1304 }
1305 nargs = a-arglist;
1306 AT.WorkPointer = a;
1307/*
1308 Now sort the list. Bubble type sort.
1309*/
1310 a1 = arglist; a2 = a1+1;
1311 while ( a2 < a ) {
1312 a3 = a2;
1313 while ( a3 > a1 && a3[-1] > a3[0] ) {
1314 EXCH(*a3,a3[-1]);
1315 a3--;
1316 }
1317 a2++;
1318 }
1319/*
1320 Now remove the elements from the list that do not occur exactly twice.
1321*/
1322 a1 = arglist; a2 = a1; a3 = a1+nargs;
1323 while ( a2 < a3 ) {
1324 if ( a2+1 == a3 ) { break; }
1325 else if ( a2+2 == a3 ) {
1326 if ( a2[0] == a2[1] ) { *a1++ = a2[0]; }
1327 break;
1328 }
1329 else {
1330 if ( a2[0] != a2[1] ) { a2++; }
1331 else if ( a2[0] != a2[2] ) {
1332 *a1++ = a2[0]; a2 += 2;
1333 }
1334 else {
1335 a2++; while ( a2 < a3 && a2[-1] == a2[0] ) a2++;
1336 }
1337 }
1338 }
1339 nargs = a1-arglist;
1340/*
1341 Now we need to redo the list of vertices and remove the elements
1342 that are not in our arglist.
1343*/
1344 do {
1345 action = 0;
1346 for ( i = -2; i < nvert; i++ ) {
1347 vv = v = AT.pWorkSpace[vert+i];
1348 v++;
1349 for ( j = 1; j < vv[0]; j++ ) {
1350 for ( jj = 0; jj < nargs; jj++ ) {
1351 if ( *v == arglist[jj] ) break;
1352 }
1353 if ( jj >= nargs ) { /* was not in the list */
1354 vv[0] = vv[0]-1;
1355 for ( jj = j; jj < vv[0]; jj++ ) vv[jj] = vv[jj+1];
1356 }
1357 v++;
1358 }
1359 }
1360/*
1361 Next we need to remove vertices that have only one object remaining.
1362 Also, the object can be removed from arglist. After that we go back
1363 to clean up the list.
1364*/
1365 for ( i = -2; i < nvert; i++ ) {
1366 vv = AT.pWorkSpace[vert+i];
1367 if ( vv[0] == 1 ) {
1368 AT.pWorkSpace[vert+i] = AT.pWorkSpace[vert+nvert-1];
1369 nvert--; i--; continue;
1370 }
1371 else if ( vv[0] == 2 && i >= 0 ) {
1372 for ( j = 0; j < nargs; j++ ) {
1373 if ( arglist[j] == vv[1] ) break;
1374 }
1375 while ( j < nargs-1 ) arglist[j] = arglist[j+1];
1376 nargs--;
1377 AT.pWorkSpace[vert+i] = AT.pWorkSpace[vert+nvert-1];
1378 nvert--; i--;
1379 action = 1;
1380 }
1381 }
1382 } while ( action );
1383/*
1384 Now we have a clean list of connections and we can start building
1385 the path. We start with summing over all objects in vert1.
1386*/
1387 path = AT.WorkPointer; npath = 0;
1388 AT.WorkPointer += nvert+8;
1389
1390 t = AT.pWorkSpace[vert1];
1391 if ( *t >= 2 ) {
1392 for ( i = 1; i < *t; i++ ) {
1393 AT.pWorkSpace[vert+nvert] = t;
1394 path[npath++] = t[i];
1395 numgenerated += GenPaths(BHEAD term,level,vert,nvert,arglist,nargs,path,npath);
1396 npath--;
1397 }
1398 }
1399 AT.WorkPointer = oldworkpointer;
1400 AT.pWorkPointer = oldpworkpointer;
1401 if ( numgenerated == 0 && option2 != 0 ) return(Generator(BHEAD term,level));
1402 return(0);
1403}
1404
1405/*
1406 #] AllPaths :
1407 #[ GenPaths :
1408
1409 We enter with an open line in path[npath-1]
1410 We traverse though the vertices until we reach vert-1, the endpoint.
1411*/
1412
1413LONG GenPaths(PHEAD WORD *term, WORD level, LONG vert, WORD nvert,
1414 WORD *arglist, WORD nargs, WORD *path, WORD npath)
1415{
1416 LONG numgenerated = 0;
1417 WORD *t, *vpartner, *v, *vv;
1418 WORD i, j;
1419/*
1420 Check whether path[npath-1] is part of the endpoint.
1421*/
1422 t = AT.pWorkSpace[vert-1];
1423 for ( i = 1; i < *t; i++ ) {
1424 if ( t[i] == path[npath-1] ) { /* Got a path! */
1425 PathOutput(BHEAD term,level,path,npath);
1426 numgenerated++;
1427 return(numgenerated);
1428 }
1429 }
1430/*
1431 Start with finding the partner.
1432*/
1433 for ( i = 0; i < nvert; i++ ) {
1434 vpartner = AT.pWorkSpace[vert+i];
1435 for ( j = 0; j < npath; j++ ) {
1436 if ( vpartner == AT.pWorkSpace[vert+nvert+j] ) break;
1437 }
1438 if ( j < npath ) continue;
1439 v = vpartner+1; vv = vpartner + *vpartner;
1440 while ( v < vv ) {
1441 if ( *v == path[npath-1] ) {
1442/*
1443 Found the partner.
1444 Now we can sum over the remaining permitted arguments of this vertex.
1445*/
1446 v = vpartner + 1;
1447 while ( v < vv ) {
1448/*
1449 *v should be in arglist.
1450*/
1451 for ( j = 0; j < nargs; j++ ) {
1452 if ( *v == arglist[j] ) break;
1453 }
1454 if ( j >= nargs ) { v++; continue; }
1455/*
1456 *v should not be in path.
1457*/
1458 for ( j = 0; j < npath; j++ ) {
1459 if ( *v == path[j] ) break;
1460 }
1461 if ( j >= npath ) {
1462 AT.pWorkSpace[vert+nvert+npath] = vpartner;
1463 path[npath++] = *v;
1464 numgenerated += GenPaths(BHEAD term,level,vert,nvert,
1465 arglist,nargs,path,npath);
1466 npath--;
1467 }
1468 v++;
1469 }
1470 return(numgenerated);
1471 }
1472 v++;
1473 }
1474 }
1475/*
1476 The partner is in a vertex we have finished before.
1477*/
1478 return(numgenerated);
1479}
1480
1481/*
1482 #] GenPaths :
1483 #[ PathOutput :
1484*/
1485
1486void PathOutput(PHEAD WORD *term, WORD level, WORD *path, WORD npath)
1487{
1488 CBUF *C = cbuf+AM.rbufnum;
1489 WORD pathfun = C->lhs[level][4]; /* The output function */
1490 WORD option1 = C->lhs[level][5]; /* type of argument */
1491 WORD *tstop, *tstop1, *t, *tt;
1492 WORD *outterm;
1493 WORD i;
1494 tstop1 = term+*term; tstop = tstop1 - ABS(tstop1[-1]);
1495 outterm = AT.WorkPointer;
1496 tt = outterm; t = term;
1497 while ( t < tstop ) *tt++ = *t++;
1498 *tt++ = pathfun;
1499 if ( functions[pathfun-FUNCTION].spec == TENSORFUNCTION ) {
1500 *tt++ = FUNHEAD+npath;
1501 FILLFUN(tt)
1502 if ( option1 == -VECTOR ) {
1503 for ( i = 0; i < npath; i++ ) *tt++ = path[i]+AM.OffsetVector;
1504 }
1505 else {
1506 for ( i = 0; i < npath; i++ ) *tt++ = path[i];
1507 }
1508 }
1509 else {
1510 *tt++ = FUNHEAD+npath*2;
1511 FILLFUN(tt)
1512 for ( i = 0; i< npath; i++ ) {
1513 *tt++ = option1;
1514 if ( option1 == -VECTOR ) *tt++ = path[i] + AM.OffsetVector;
1515 else *tt++ = path[i];
1516 }
1517 }
1518 while ( t < tstop1 ) *tt++ = *t++;
1519 *outterm = tt - outterm;
1520 AT.WorkPointer = tt;
1521 if ( Generator(BHEAD outterm,level) ) {
1522 MesCall("PathOutput");
1523 Terminate(-1);
1524 }
1525 AT.WorkPointer = outterm;
1526}
1527
1528
1529
1530/*
1531 #] PathOutput :
1532 #[ AllOnePI :
1533
1534 We assume a graph that has loops and is OnePI.
1535 This routine initializes the recursion that creates all onePI subgraphs.
1536
1537 Algorithm:
1538 a: have a routine that removes all bridges: RemoveBridges
1539 b: output an empty diagram.
1540 c: output the diagram itself
1541 d: cut a line, followed by removing all bridges.
1542 e: if the diagram still has lines, go to c, else go to the next line in d.
1543 In order to avoid factorial blowup, we need to prune the tree as fast as possible.
1544
1545 1: list of lines to be cut.
1546 2: once we have tried a line and removed its bridges, we do not have to
1547 try this line deeper in the tree, neither its bridges.
1548
1549
1550 1 2 3
1551 cut 1. x
1552 cut 2 -> bridge 3 x
1553 cut 2. x
1554 cut 3 -> bridge 1 ???? mag niet
1555 cut 3. x
1556 Hence 1,2,3,4,5,6,7
1557 1 still to go 2,3,4,5,6,7
1558 2 still to go 3,4,5,6,7
1559 3 still to go 4,5,6,7
1560 etc.
1561 bridges: if x is bridge, we take x from the list_to_go.
1562 If we have a bridge that is a number lower than the list_to_go we skip this possibility.
1563 We reach the end when the list_to_go is empty.
1564*/
1565
1566WORD AllOnePI(WORD *term,WORD level)
1567{
1568 CBUF *C = cbuf+AM.rbufnum;
1569 WORD vcode = C->lhs[level][2]; /* The vertex function */
1570 WORD option1 = C->lhs[level][4]; /* type of argument */
1571/*
1572 First we have to collect all relevant information about the diagram.
1573 We should start with removing bridges to make the real starting point onePI.
1574*/
1575 DUMMYUSE(term)
1576 DUMMYUSE(vcode)
1577 DUMMYUSE(option1)
1578 return(0);
1579}
1580
1581/*
1582 #] AllOnePI :
1583 #[ RemoveBridges :
1584*/
1585
1586int RemoveBridges(void)
1587{
1588 return(0);
1589}
1590
1591/*
1592 #] RemoveBridges :
1593 #[ TakeOneLine :
1594*/
1595
1596int TakeOneLine(WORD*term,WORD level)
1597{
1598 DUMMYUSE(term)
1599 DUMMYUSE(level)
1600 return(0);
1601}
1602
1603/*
1604 #] TakeOneLine :
1605 #[ OutputOnePI :
1606*/
1607
1608int OutputOnePI(PHEAD WORD *term,WORD level)
1609{
1610 return(Generator(BHEAD term,level));
1611}
1612
1613/*
1614 #] OutputOnePI :
1615*/
1616
Generator
int Generator(PHEAD WORD *, WORD)
Definition proces.c:3249
CbUf::lhs
WORD ** lhs
Definition structs.h:974
Generated by doxygen 1.9.8