FORM 4.3
lus.c
Go to the documentation of this file.
1
7/* #[ License : */
8/*
9 * Copyright (C) 1984-2022 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*/