57int Lus(WORD *term, WORD funnum, WORD loopsize, WORD numargs, WORD outfun, WORD mode)
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);
69 if ( ( ten = functions[funnum-FUNCTION].spec ) >= TENSORFUNCTION ) {
71 if ( *t == funnum && t[1] == FUNHEAD+numargs ) { nfun++; }
78 i = 0; m = t+FUNHEAD; t += t[1];
79 while ( m < t ) { i++; NEXTARG(m) }
80 if ( i == numargs ) nfun++;
85 if ( loopsize < 0 ) minloopsize = 2;
86 else minloopsize = loopsize;
87 if ( funnum < minloopsize )
return(0);
88 if ( ((functions[funnum-FUNCTION].symmetric) & ~REVERSEORDER) == ANTISYMMETRIC ) sign = 1;
89 if ( mode == 1 || mode < 0 ) {
90 ten2 = functions[outfun-FUNCTION].spec >= TENSORFUNCTION;
96 if ( AN.numflocs < funnum ) {
97 if ( AN.funlocs ) M_free(AN.funlocs,
"Lus: AN.funlocs");
99 AN.funlocs = (WORD **)Malloc1(
sizeof(WORD *)*AN.numflocs,
"Lus: AN.funlocs");
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");
109 alist = AN.funargs; loc = AN.funlocs;
111 if ( ten >= TENSORFUNCTION ) {
112 while ( t < tstop ) {
113 if ( *t == funnum && t[1] == FUNHEAD+numargs ) {
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--;
124 while ( --j >= 0 ) *alist++ = -1;
131 while ( t < tstop ) {
132 if ( *t == funnum ) {
134 i = 0; m = t+FUNHEAD; t += t[1];
135 while ( m < t ) { i++; NEXTARG(m) }
136 if ( i == numargs ) {
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--;
144 else if ( ten2 >= TENSORFUNCTION && *m != -INDEX
145 && *m != -VECTOR && *m != -MINVECTOR &&
146 ( *m != -SNUMBER || *m < 0 || *m >= AM.OffsetIndex ) ) {
154 while ( --i >= 0 ) *alist++ = -1;
160 if ( nfun < minloopsize )
return(0);
169 alist = AN.funargs; totnum = numargs*nfun;
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");
176 i = totnum; n = 0; wi = AN.funinds;
178 if ( *alist >= 0 ) { n++; *wi++ = *alist; *wi++ = 1; }
181 n = SortTheList(AN.funinds,n);
184 for ( i = 0; i < nfun; i++ ) {
185 alist = AN.funargs + i*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;
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++; }
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;
211 (AN.funinds[2*medi+1])--;
212 if ( AN.funinds[2*medi+1] == 1 ) action++;
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++;
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;
228 if ( j >= totnum ) alist[i] = -1;
232 for ( i = 0; i < nfun; i++ ) {
233 alist = AN.funargs + i*numargs;
235 for ( k = 0; k < n; k++ ) {
236 if ( alist[k] < 0 ) { alist[k--] = alist[--n]; alist[n] = -1; }
239 if ( n == 1 ) { j = alist[0]; }
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];
246 for ( k = 0, jj = 0, wi = AN.funargs; k < totnum; k++, wi++ ) {
247 if ( *wi == j ) { jj++;
if ( jj > 1 )
break; }
250 for ( k = 0, wi = AN.funargs; k < totnum; k++, wi++ ) {
251 if ( *wi == j ) { *wi = -1; action = 1; }
259 if ( nfun < minloopsize )
return(0);
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;
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;
276 else if ( FindLus(0,loopsize,AN.tohunt) ) { L = loopsize;
goto Success; }
279 AN.nargs = numargs; AN.numoffuns = nfun;
280 if ( loopsize < 0 ) {
282 if ( loopsize < -1 ) { k = -loopsize-1;
if ( k > nfun ) k = nfun; }
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];
305 if ( mode == 0 || mode > 1 )
return(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;
315 while ( t < tstop ) {
316 if ( *t == -1 )
break;
319 while ( m < t ) *w++ = *m++;
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] ) {
333 jj = m[k]; m[k] = m[0]; m[0] = jj;
339 for ( k = 1; k < numargs; k++ ) {
340 if ( m[k] == alist[1] ) {
342 jj = m[k]; m[k] = m[1]; m[1] = jj;
349 for ( k = 2; k < numargs; k++ ) *w++ = *m++;
354 for ( i = 0; i < L; i++ ) {
355 alist = wi + i*numargs;
358 for ( k = 0; k < numargs; k++ ) {
359 if ( *m == -INDEX && m[1] == alist[0] ) {
361 if ( ( k & 1 ) != 0 ) sign = -sign;
365 t2 = m+2; t1 = m; t3 = tt+FUNHEAD;
366 while ( t1 > t3 ) { *--t2 = *--t1; }
367 t3[0] = -INDEX; t3[1] = alist[0];
373 m = tt + FUNHEAD + 2;
374 for ( k = 1; k < numargs; k++ ) {
375 if ( *m == -INDEX && m[1] == alist[1] ) {
377 if ( ( k & 1 ) == 0 ) sign = -sign;
381 t2 = m+2; t1 = m; t3 = tt+FUNHEAD+2;
382 while ( t1 > t3 ) { *--t2 = *--t1; }
383 t3[0] = -INDEX; t3[1] = alist[1];
393 t1 = tt + FUNHEAD + 4;
396 if ( *t1 == -INDEX || *t1 == -VECTOR ) {
397 *w++ = t1[1]; t1 += 2;
399 else if ( *t1 == -MINVECTOR ) {
400 *w++ = t1[1]; t1 += 2; sign2 = -sign2;
402 else if ( ( *t1 == -SNUMBER ) && ( t1[1] >= 0 ) && ( t1[1] < AM.OffsetIndex ) ) {
403 *w++ = t1[1]; t1 += 2; sign2 = -sign2;
406 MLOCK(ErrorMessageLock);
407 MesPrint(
"Illegal attempt to use a non-index-like argument in a tensor in ReplaceLoop statement");
408 MUNLOCK(ErrorMessageLock);
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] ) {
423 jj = m[k]; m[k] = m[0]; m[0] = jj;
429 for ( k = 1; k < numargs; k++ ) {
430 if ( m[k] == alist[1] ) {
432 jj = m[k]; m[k] = m[1]; m[1] = jj;
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; }
449 for ( i = 0; i < L; i++ ) {
450 alist = wi + i*numargs;
453 for ( k = 0; k < numargs; k++ ) {
454 if ( *m == -INDEX && m[1] == alist[0] ) {
456 if ( ( k & 1 ) != 0 ) sign = -sign;
460 t2 = m+2; t1 = m; t3 = tt+FUNHEAD;
461 while ( t1 > t3 ) { *--t2 = *--t1; }
462 t3[0] = -INDEX; t3[1] = alist[0];
468 m = tt + FUNHEAD + 2;
469 for ( k = 1; k < numargs; k++ ) {
470 if ( *m == -INDEX && m[1] == alist[1] ) {
472 if ( ( k & 1 ) == 0 ) sign = -sign;
476 t2 = m+2; t1 = m; t3 = tt+FUNHEAD+2;
477 while ( t1 > t3 ) { *--t2 = *--t1; }
478 t3[0] = -INDEX; t3[1] = alist[1];
487 t1 = tt + FUNHEAD + 4;
489 while ( t1 < t2 ) *w++ = *t1++;
494 while ( t < tstop ) {
495 if ( *t == -1 ) { t += t[1];
continue; }
499 tstop = term + *term;
500 while ( t < tstop ) *w++ = *t++;
501 if ( sign < 0 ) w[-1] = -w[-1];
502 i = w - AT.WorkPointer;
504 t = term; w = AT.WorkPointer;
515int FindLus(
int from,
int level,
int openindex)
518 int i, j, k, jj, *alist, *blist, *w, *m, partner;
519 WORD **loc = AN.funlocs, *wor;
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 ) {
529 alist = AN.funargs + from*AN.nargs;
530 alist[0] = openindex; alist[1] = AN.tohunt;
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 ) {
544 wor = loc[i]; loc[i] = loc[from]; loc[from] = wor;
545 blist = w = AN.funargs + from*AN.nargs;
548 while ( --k >= 0 ) { jj = *m; *m++ = *w; *w++ = jj; }
551 for ( k = 0; k < AN.nargs; k++ ) {
552 if ( k == j || blist[k] < 0 )
continue;
554 if ( FindLus(from+1,level-1,partner) ) {
555 blist[0] = openindex; blist[1] = partner;
560 wor = loc[i]; loc[i] = loc[from]; loc[from] = wor;
561 w = AN.funargs + from*AN.nargs;
564 while ( --k >= 0 ) { jj = *m; *m++ = *w; *w++ = jj; }
578int SortTheList(
int *slist,
int num)
581 int i, nleft, nright, *t1, *t2, *t3, *rlist;
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;
590 slist[1] += slist[3];
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");
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 ) {
607 *t3++ = *t1++; *t3++ = *t1++; nleft--;
609 else if ( *t1 > *t2 ) {
610 *t3++ = *t2++; *t3++ = *t2++; nright--;
613 *t3++ = *t1++; t2++; *t3++ = (*t1++) + (*t2++); i--;
617 while ( --nleft >= 0 ) { *t3++ = *t1++; *t3++ = *t1++; }
618 while ( --nright >= 0 ) { *t3++ = *t2++; *t3++ = *t2++; }