FORM  4.2.1
lus.c
Go to the documentation of this file.
1 
7 /* #[ License : */
8 /*
9  * Copyright (C) 1984-2017 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 
57 int 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);
304 Success:;
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 
515 int 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 
578 int 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 */