FORM 4.3
reshuf.c
Go to the documentation of this file.
1
9/* #[ License : */
10/*
11 * Copyright (C) 1984-2022 J.A.M. Vermaseren
12 * When using this file you are requested to refer to the publication
13 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
14 * This is considered a matter of courtesy as the development was paid
15 * for by FOM the Dutch physics granting agency and we would like to
16 * be able to track its scientific use to convince FOM of its value
17 * for the community.
18 *
19 * This file is part of FORM.
20 *
21 * FORM is free software: you can redistribute it and/or modify it under the
22 * terms of the GNU General Public License as published by the Free Software
23 * Foundation, either version 3 of the License, or (at your option) any later
24 * version.
25 *
26 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
27 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
28 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
29 * details.
30 *
31 * You should have received a copy of the GNU General Public License along
32 * with FORM. If not, see <http://www.gnu.org/licenses/>.
33 */
34/* #] License : */
35#define NEWCODE
36/*
37 #[ Includes : reshuf.c
38*/
39
40#include "form3.h"
41
42/*
43 #] Includes :
44 #[ Reshuf :
45
46 Routines to rearrange dummy indices, so that
47 a: The notation becomes reasonably unique (the perfect job
48 may consume very much time).
49 b: The value of AR.CurDum is reset.
50
51 Also some routines are needed to aid in the reading of stored
52 expressions. Also in those expressions there can be dummy
53 indices, and there should be no conflict with the already
54 existing dummies.
55
56 #[ ReNumber :
57
58 Reads the term, tests for dummies, and puts them in order.
59 Note that this is kind of a first order approximation.
60 There is quite some room to make this routine 'smart'
61 First order:
62 First index will be lowest, second will be next etc.
63 Second order:
64 Functions with more than one index and symmetry properties
65 have some look ahead to see which index is the first to
66 find its partner.
67 Third order:
68 Take the ordering of the functions into account.
69 Fourth order:
70 Try all permutations and see which one gives the 'minimal' term.
71 Currently we use only the first order.
72
73 We need a scratch array for the numbers we find, and one for
74 the addresses at which these numbers are.
75 We can use the space for the Scrat arrays. There are 13 of those
76 and each is AM.MaxTal UWORDs long.
77
78*/
79
80WORD ReNumber(PHEAD WORD *term)
81{
82 GETBIDENTITY
83 WORD *d, *e, **p, **f;
84 WORD n, i, j, old;
85 AN.DumFound = AN.RenumScratch;
86 AN.DumPlace = AN.PoinScratch;
87 AN.DumFunPlace = AN.FunScratch;
88 AN.NumFound = 0;
89 FunLevel(BHEAD term);
90 d = AN.RenumScratch;
91 p = AN.PoinScratch;
92 f = AN.FunScratch;
93/*
94 Now the first level sorting.
95*/
96 i = AN.IndDum;
97 n = AN.NumFound;
98 while ( --n >= 0 ) {
99 if ( *d > 0 ) {
100 old = **p;
101 **p = ++i;
102 if ( *f ) **f |= DIRTYSYMFLAG;
103 e = d;
104 e++;
105 for ( j = 1; j <= n; j++ ) {
106 if ( *e && *(p[j]) == old ) {
107 *(p[j]) = i;
108 if ( f[j] ) *(f[j]) |= DIRTYSYMFLAG;
109 *e = 0;
110 }
111 e++;
112 }
113 }
114 p++;
115 d++;
116 f++;
117 }
118 return(i);
119}
120
121/*
122 #] ReNumber :
123 #[ FunLevel :
124
125 Does one term in determining where the dummies are.
126 Made to work recursively for functions.
127
128*/
129
130VOID FunLevel(PHEAD WORD *term)
131{
132 GETBIDENTITY
133 WORD *t, *tstop, *r, *fun;
134 WORD *m;
135 t = r = term;
136 r += *r;
137 tstop = r - ABS(r[-1]);
138 t++;
139 if ( t < tstop ) do {
140 r = t + t[1];
141 switch ( *t ) {
142 case SYMBOL:
143 case DOTPRODUCT:
144 break;
145 case VECTOR:
146 t += 3;
147 do {
148 if ( *t > AN.IndDum ) {
149 if ( AN.NumFound >= AN.MaxRenumScratch ) AdjustRenumScratch(BHEAD0);
150 AN.NumFound++;
151 *AN.DumFound++ = *t;
152 *AN.DumPlace++ = t;
153 *AN.DumFunPlace++ = 0;
154 }
155 t += 2;
156 } while ( t < r );
157 break;
158 case SUBEXPRESSION:
159/*
160 Still must hunt down the wildcards(?)
161*/
162 break;
163 case GAMMA:
164 t += FUNHEAD-2;
165 /* fall through */
166 case DELTA:
167 case INDEX:
168 t += 2;
169 while ( t < r ) {
170 if ( *t > AN.IndDum ) {
171 if ( AN.NumFound >= AN.MaxRenumScratch ) AdjustRenumScratch(BHEAD0);
172 AN.NumFound++;
173 *AN.DumFound++ = *t;
174 *AN.DumPlace++ = t;
175 *AN.DumFunPlace++ = 0;
176 }
177 t++;
178 }
179 break;
180 case HAAKJE:
181 case EXPRESSION:
182 case SNUMBER:
183 case LNUMBER:
184 break;
185 default:
186 if ( *t < FUNCTION ) {
187 MLOCK(ErrorMessageLock);
188 MesPrint("Unexpected code in ReNumber");
189 MUNLOCK(ErrorMessageLock);
190 Terminate(-1);
191 }
192 fun = t+2;
193 if ( *t >= FUNCTION && functions[*t-FUNCTION].spec
194 >= TENSORFUNCTION ) {
195 t += FUNHEAD;
196 while ( t < r ) {
197 if ( *t > AN.IndDum ) {
198 if ( AN.NumFound >= AN.MaxRenumScratch ) AdjustRenumScratch(BHEAD0);
199 AN.NumFound++;
200 *AN.DumFound++ = *t;
201 *AN.DumPlace++ = t;
202 *AN.DumFunPlace++ = fun;
203 }
204 t++;
205 }
206 break;
207 }
208
209 t += FUNHEAD;
210 while ( t < r ) {
211 if ( *t > 0 ) {
212
213 /* General function. Enter 'recursion'. */
214
215 m = t + *t;
216 t += ARGHEAD;
217 while ( t < m ) {
218 FunLevel(BHEAD t);
219 t += *t;
220 }
221 }
222 else {
223 if ( *t == -INDEX ) {
224 t++;
225 if ( *t >= AN.IndDum ) {
226 if ( AN.NumFound >= AN.MaxRenumScratch ) AdjustRenumScratch(BHEAD0);
227 AN.NumFound++;
228 *AN.DumFound++ = *t;
229 *AN.DumPlace++ = t;
230 *AN.DumFunPlace++ = fun;
231 }
232 t++;
233 }
234 else if ( *t <= -FUNCTION ) t++;
235 else t += 2;
236 }
237 }
238 break;
239 }
240 t = r;
241 } while ( t < tstop );
242}
243
244/*
245 #] FunLevel :
246 #[ DetCurDum :
247
248 We look for indices in the range AM.IndDum to AM.IndDum+MAXDUMMIES.
249 The maximum value is returned.
250*/
251
252WORD DetCurDum(PHEAD WORD *t)
253{
254 GETBIDENTITY
255 WORD maxval = AN.IndDum;
256 WORD maxtop = AM.IndDum + WILDOFFSET;
257 WORD *tstop, *m, *r, i;
258 tstop = t + *t - 1;
259 tstop -= ABS(*tstop);
260 t++;
261 while ( t < tstop ) {
262 if ( *t == VECTOR ) {
263 m = t + 3;
264 t += t[1];
265 while ( m < t ) {
266 if ( *m > maxval && *m < maxtop ) maxval = *m;
267 m += 2;
268 }
269 }
270 else if ( *t == DELTA || *t == INDEX ) {
271 m = t + 2;
272Singles:
273 t += t[1];
274 while ( m < t ) {
275 if ( *m > maxval && *m < maxtop ) maxval = *m;
276 m++;
277 }
278 }
279 else if ( *t >= FUNCTION ) {
280 if ( functions[*t-FUNCTION].spec >= TENSORFUNCTION ) {
281 m = t + FUNHEAD;
282 goto Singles;
283 }
284 r = t + FUNHEAD;
285 t += t[1];
286 while ( r < t ) { /* The arguments */
287 if ( *r < 0 ) {
288 if ( *r <= -FUNCTION ) r++;
289 else if ( *r == -INDEX ) {
290 if ( r[1] > maxval && r[1] < maxtop ) maxval = r[1];
291 r += 2;
292 }
293 else r += 2;
294 }
295 else {
296 m = r + ARGHEAD;
297 r += *r;
298 while ( m < r ) { /* Terms in the argument */
299 i = DetCurDum(BHEAD m);
300 if ( i > maxval && i < maxtop ) maxval = i;
301 m += *m;
302 }
303 }
304 }
305 }
306 else {
307 t += t[1];
308 }
309 }
310 return(maxval);
311}
312
313/*
314 #] DetCurDum :
315 #[ FullRenumber :
316
317 Does a full renumbering. May be slow if there are many indices.
318 par = 1 Goes with a factorial!
319 par = 0 All single exchanges only till there is no more improvement.
320 Notice that there is a hole in the defense with respect to
321 arguments inside functions inside functions.
322*/
323
324int FullRenumber(PHEAD WORD *term, WORD par)
325{
326 GETBIDENTITY
327 WORD *d, **p, **f, *w, *t, *best, *stac, *perm, a, *termtry;
328 WORD n, i, j, k, ii;
329 WORD *oldworkpointer = AT.WorkPointer;
330 n = ReNumber(BHEAD term) - AM.IndDum;
331 if ( n <= 1 ) return(0);
332 Normalize(BHEAD term);
333 if ( *term == 0 ) return(0);
334 n = ReNumber(BHEAD term) - AM.IndDum;
335 d = AN.RenumScratch;
336 p = AN.PoinScratch;
337 f = AN.FunScratch;
338 if ( AT.WorkPointer < term + *term ) AT.WorkPointer = term + *term;
339 k = AN.NumFound;
340 best = w = AT.WorkPointer; t = term;
341 for ( i = *term; i > 0; i-- ) *w++ = *t++;
342 AT.WorkPointer = w;
343 Normalize(BHEAD best);
344 AT.WorkPointer = w = best + *best;
345 stac = w+100;
346 perm = stac + n + 1;
347 termtry = perm + n + 1;
348 for ( i = 1; i <= n; i++ ) perm[i] = i + AM.IndDum;
349 for ( i = 1; i <= n; i++ ) stac[i] = i;
350 for ( i = 0; i < k; i++ ) d[i] = *(p[i]) - AM.IndDum;
351 if ( par == 0 ) { /* All single exchanges */
352 for ( i = 1; i < n; i++ ) {
353 for ( j = i+1; j <= n; j++ ) {
354 a = perm[j]; perm[j] = perm[i]; perm[i] = a;
355 for ( ii = 0; ii < k; ii++ ) {
356 *(p[ii]) = perm[d[ii]];
357 if ( f[ii] ) *(f[ii]) |= DIRTYSYMFLAG;
358 }
359 t = term; w = termtry;
360 for ( ii = 0; ii < *term; ii++ ) *w++ = *t++;
361 AT.WorkPointer = w;
362 if ( Normalize(BHEAD termtry) == 0 ) {
363 if ( *termtry == 0 ) goto Return0;
364 if ( ( ii = CompareTerms(termtry,best,0) ) > 0 ) {
365 t = termtry; w = best;
366 for ( ii = 0; ii < *termtry; ii++ ) *w++ = *t++;
367 i = 0; break; /* restart from beginning */
368 }
369 else if ( ii == 0 && CompCoef(termtry,best) != 0 )
370 goto Return0;
371 }
372 /* if no success, set back */
373 a = perm[j]; perm[j] = perm[i]; perm[i] = a;
374 }
375 }
376 }
377 else if ( par == 1 ) { /* all permutations */
378 j = n-1;
379 for(;;) {
380 if ( stac[j] == n ) {
381 a = perm[j]; perm[j] = perm[n]; perm[n] = a;
382 stac[j] = j;
383 j--;
384 if ( j <= 0 ) break;
385 continue;
386 }
387 if ( j != stac[j] ) {
388 a = perm[j]; perm[j] = perm[stac[j]]; perm[stac[j]] = a;
389 }
390 (stac[j])++;
391 a = perm[j]; perm[j] = perm[stac[j]]; perm[stac[j]] = a;
392 {
393 for ( i = 0; i < k; i++ ) {
394 *(p[i]) = perm[d[i]];
395 if ( f[i] ) *(f[i]) |= DIRTYSYMFLAG;
396 }
397 t = term; w = termtry;
398 for ( i = 0; i < *term; i++ ) *w++ = *t++;
399 AT.WorkPointer = w;
400 if ( Normalize(BHEAD termtry) == 0 ) {
401 if ( *termtry == 0 ) goto Return0;
402 if ( ( ii = CompareTerms(termtry,best,0) ) > 0 ) {
403 t = termtry; w = best;
404 for ( i = 0; i < *termtry; i++ ) *w++ = *t++;
405 }
406 else if ( ii == 0 && CompCoef(termtry,best) != 0 )
407 goto Return0;
408 }
409 }
410 if ( j < n-1 ) { j = n-1; }
411 }
412 }
413 t = term; w = best;
414 n = *best;
415 for ( i = 0; i < n; i++ ) *t++ = *w++;
416 AT.WorkPointer = oldworkpointer;
417 return(0);
418Return0:
419 *term = 0;
420 AT.WorkPointer = oldworkpointer;
421 return(0);
422}
423
424/*
425 #] FullRenumber :
426 #[ MoveDummies :
427
428 Routine shifts the dummy indices by an amount 'shift'.
429 It is an adaptation of DetCurDum.
430 Needed for = ...*expression^power*...
431 in which expression contains dummy indices.
432 Note that this code should have been in ver1 already and has
433 always been missing. Routine made 29-jan-2007.
434*/
435
436VOID MoveDummies(PHEAD WORD *term, WORD shift)
437{
438 GETBIDENTITY
439 WORD maxval = AN.IndDum;
440 WORD maxtop = AM.IndDum + WILDOFFSET;
441 WORD *tstop, *m, *r;
442 tstop = term + *term - 1;
443 tstop -= ABS(*tstop);
444 term++;
445 while ( term < tstop ) {
446 if ( *term == VECTOR ) {
447 m = term + 3;
448 term += term[1];
449 while ( m < term ) {
450 if ( *m > maxval && *m < maxtop ) *m += shift;
451 m += 2;
452 }
453 }
454 else if ( *term == DELTA || *term == INDEX ) {
455 m = term + 2;
456Singles:
457 term += term[1];
458 while ( m < term ) {
459 if ( *m > maxval && *m < maxtop ) *m += shift;
460 m++;
461 }
462 }
463 else if ( *term >= FUNCTION ) {
464 if ( functions[*term-FUNCTION].spec >= TENSORFUNCTION ) {
465 m = term + FUNHEAD;
466 goto Singles;
467 }
468 r = term + FUNHEAD;
469 term += term[1];
470 while ( r < term ) { /* The arguments */
471 if ( *r < 0 ) {
472 if ( *r <= -FUNCTION ) r++;
473 else if ( *r == -INDEX ) {
474 if ( r[1] > maxval && r[1] < maxtop ) r[1] += shift;
475 r += 2;
476 }
477 else r += 2;
478 }
479 else {
480 m = r + ARGHEAD;
481 r += *r;
482 while ( m < r ) { /* Terms in the argument */
483 MoveDummies(BHEAD m,shift);
484 m += *m;
485 }
486 }
487 }
488 }
489 else {
490 term += term[1];
491 }
492 }
493}
494
495/*
496 #] MoveDummies :
497 #[ AdjustRenumScratch :
498
499 Extends the buffer for number of dummies that can be found in
500 a term. Originally we had a fixed buffer at size 300 in the AN
501 struct. Thomas Hahn ran out of that. Hence we have now made it
502 a dynamical buffer.
503 Note that the pointers used in FunLevel need adjustment as well.
504*/
505
506void AdjustRenumScratch(PHEAD0)
507{
508 GETBIDENTITY
509 WORD newsize;
510 int i;
511 WORD **newpoin, *newnum;
512 if ( AN.MaxRenumScratch == 0 ) newsize = 100;
513 else newsize = AN.MaxRenumScratch*2;
514 if ( newsize > MAXPOSITIVE/2 ) newsize = MAXPOSITIVE/2+1;
515
516 newpoin = (WORD **)Malloc1(newsize*sizeof(WORD *),"PoinScratch");
517 for ( i = 0; i < AN.NumFound; i++ ) newpoin[i] = AN.PoinScratch[i];
518 for ( ; i < newsize; i++ ) newpoin[i] = 0;
519 if ( AN.PoinScratch ) M_free(AN.PoinScratch,"PoinScratch");
520 AN.PoinScratch = newpoin;
521 AN.DumPlace = newpoin + AN.NumFound;
522
523 newpoin = (WORD **)Malloc1(newsize*sizeof(WORD *),"FunScratch");
524 for ( i = 0; i < AN.NumFound; i++ ) newpoin[i] = AN.FunScratch[i];
525 for ( ; i < newsize; i++ ) newpoin[i] = 0;
526 if ( AN.FunScratch ) M_free(AN.FunScratch,"FunScratch");
527 AN.FunScratch = newpoin;
528 AN.DumFunPlace = newpoin + AN.NumFound;
529
530 newnum = (WORD *)Malloc1(newsize*sizeof(WORD),"RenumScratch");
531 for ( i = 0; i < AN.NumFound; i++ ) newnum[i] = AN.RenumScratch[i];
532 for ( ; i < newsize; i++ ) newnum[i] = 0;
533 if ( AN.RenumScratch ) M_free(AN.RenumScratch,"RenumScratch");
534 AN.RenumScratch = newnum;
535 AN.DumFound = newnum + AN.NumFound;
536
537 AN.MaxRenumScratch = newsize;
538}
539
540/*
541 #] AdjustRenumScratch :
542 #] Reshuf :
543 #[ Count :
544 #[ CountDo :
545
546 This function executes the counting action in a count
547 operation. The return value is the count of the term.
548 Input is the term and a pointer to the instruction.
549
550*/
551
552WORD CountDo(WORD *term, WORD *instruct)
553{
554 WORD *m, *r, i, j, count = 0;
555 WORD *stopper, *tstop, *r1 = 0, *r2 = 0;
556 m = instruct;
557 stopper = m + m[1];
558 instruct += 3;
559 tstop = term + *term; tstop -= ABS(tstop[-1]); term++;
560 while ( term < tstop ) {
561 switch ( *term ) {
562 case SYMBOL:
563 i = term[1] - 2;
564 term += 2;
565 while ( i > 0 ) {
566 m = instruct;
567 while ( m < stopper ) {
568 if ( *m == SYMBOL && m[2] == *term ) {
569 count += m[3] * term[1];
570 }
571 m += m[1];
572 }
573 term += 2;
574 i -= 2;
575 }
576 break;
577 case DOTPRODUCT:
578 i = term[1] - 2;
579 term += 2;
580 while ( i > 0 ) {
581 m = instruct;
582 while ( m < stopper ) {
583 if ( *m == DOTPRODUCT && (( m[2] == *term &&
584 m[3] == term[1]) || ( m[2] == term[1] &&
585 m[3] == *term )) ) {
586 count += m[4] * term[2];
587 break;
588 }
589 m += m[1];
590 }
591 m = instruct;
592 while ( m < stopper ) {
593 if ( *m == VECTOR && m[2] == *term &&
594 ( m[3] & DOTPBIT ) != 0 ) {
595 count += m[m[1]-1] * term[2];
596 }
597 m += m[1];
598 }
599 m = instruct;
600 while ( m < stopper ) {
601 if ( *m == VECTOR && m[2] == term[1] &&
602 ( m[3] & DOTPBIT ) != 0 ) {
603 count += m[m[1]-1] * term[2];
604 }
605 m += m[1];
606 }
607 term += 3;
608 i -= 3;
609 }
610 break;
611 case INDEX:
612 j = 1;
613 goto VectInd;
614 case VECTOR:
615 j = 2;
616VectInd: i = term[1] - 2;
617 term += 2;
618 while ( i > 0 ) {
619 m = instruct;
620 while ( m < stopper ) {
621 if ( *m == VECTOR && m[2] == *term &&
622 ( m[3] & VECTBIT ) != 0 ) {
623 count += m[m[1]-1];
624 }
625 m += m[1];
626 }
627 term += j;
628 i -= j;
629 }
630 break;
631 default:
632 if ( *term >= FUNCTION ) {
633 i = *term;
634 m = instruct;
635 while ( m < stopper ) {
636 if ( *m == FUNCTION && m[2] == i ) count += m[3];
637 m += m[1];
638 }
639 if ( functions[i-FUNCTION].spec >= TENSORFUNCTION ) {
640 i = term[1] - FUNHEAD;
641 term += FUNHEAD;
642 while ( i > 0 ) {
643 if ( *term < 0 ) {
644 m = instruct;
645 while ( m < stopper ) {
646 if ( *m == VECTOR && m[2] == *term &&
647 ( m[3] & FUNBIT ) != 0 ) {
648 count += m[m[1]-1];
649 }
650 m += m[1];
651 }
652 }
653 term++;
654 i--;
655 }
656 }
657 else {
658 r = term + term[1];
659 term += FUNHEAD;
660 while ( term < r ) {
661 if ( ( *term == -INDEX || *term == -VECTOR
662 || *term == -MINVECTOR ) && term[1] < MINSPEC ) {
663 m = instruct;
664 while ( m < stopper ) {
665 if ( *m == VECTOR && term[1] == m[2]
666 && ( m[3] & SETBIT ) != 0 ) {
667 r1 = SetElements + Sets[m[4]].first;
668 r2 = SetElements + Sets[m[4]].last;
669 while ( r1 < r2 ) {
670 if ( *r1 == i ) {
671 count += m[m[1]-1];
672 goto NextFF;
673 }
674 r1++;
675 }
676 }
677 m += m[1];
678 }
679NextFF:
680 term += 2;
681 }
682 else { NEXTARG(term) }
683 }
684 }
685 break;
686 }
687 else {
688 term += term[1];
689 }
690 break;
691 }
692 }
693 return(count);
694}
695
696/*
697 #] CountDo :
698 #[ CountFun :
699
700 This is the count function.
701 The return value is the count of the term.
702 Input is the term and a pointer to the count function.
703
704*/
705
706WORD CountFun(WORD *term, WORD *countfun)
707{
708 WORD *m, *r, i, j, count = 0, *instruct, *stopper, *tstop;
709 m = countfun;
710 stopper = m + m[1];
711 instruct = countfun + FUNHEAD;
712 tstop = term + *term; tstop -= ABS(tstop[-1]); term++;
713 while ( term < tstop ) {
714 switch ( *term ) {
715 case SYMBOL:
716 i = term[1] - 2;
717 term += 2;
718 while ( i > 0 ) {
719 m = instruct;
720 while ( m < stopper ) {
721 if ( *m == -SNUMBER ) { NEXTARG(m) continue; }
722 if ( *m == -SYMBOL && m[1] == *term
723 && m[2] == -SNUMBER && ( m + 2 ) < stopper ) {
724 count += m[3] * term[1]; m += 4;
725 }
726 else { NEXTARG(m) }
727 }
728 term += 2;
729 i -= 2;
730 }
731 break;
732 case DOTPRODUCT:
733 i = term[1] - 2;
734 term += 2;
735 while ( i > 0 ) {
736 m = instruct;
737 while ( m < stopper ) {
738 if ( *m == -SNUMBER ) { NEXTARG(m) continue; }
739 if ( *m == 9+ARGHEAD && m[ARGHEAD] == 9
740 && m[ARGHEAD+1] == DOTPRODUCT
741 && m[ARGHEAD+9] == -SNUMBER && ( m + ARGHEAD+9 ) < stopper
742 && (( m[ARGHEAD+3] == *term &&
743 m[ARGHEAD+4] == term[1]) ||
744 ( m[ARGHEAD+3] == term[1] &&
745 m[ARGHEAD+4] == *term )) ) {
746 count += m[ARGHEAD+10] * term[2];
747 m += ARGHEAD+11;
748 }
749 else { NEXTARG(m) }
750 }
751 m = instruct;
752 while ( m < stopper ) {
753 if ( *m == -SNUMBER ) { NEXTARG(m) continue; }
754 if ( ( *m == -VECTOR || *m == -MINVECTOR )
755 && m[1] == *term &&
756 m[2] == -SNUMBER && ( m+2 ) < stopper ) {
757 count += m[3] * term[2]; m += 4;
758 }
759 NEXTARG(m)
760 }
761 m = instruct;
762 while ( m < stopper ) {
763 if ( *m == -SNUMBER ) { NEXTARG(m) continue; }
764 if ( ( *m == -VECTOR || *m == -MINVECTOR )
765 && m[1] == term[1] &&
766 m[2] == -SNUMBER && ( m+2 ) < stopper ) {
767 count += m[3] * term[2];
768 m += 4;
769 }
770 NEXTARG(m)
771 }
772 term += 3;
773 i -= 3;
774 }
775 break;
776 case INDEX:
777 j = 1;
778 goto VectInd;
779 case VECTOR:
780 j = 2;
781VectInd: i = term[1] - 2;
782 term += 2;
783 while ( i > 0 ) {
784 m = instruct;
785 while ( m < stopper ) {
786 if ( *m == -SNUMBER ) { NEXTARG(m) continue; }
787 if ( ( *m == -VECTOR || *m == -MINVECTOR )
788 && m[1] == *term &&
789 m[2] == -SNUMBER && (m+2) < stopper ) {
790 count += m[3]; m += 4;
791 }
792 NEXTARG(m)
793 }
794 term += j;
795 i -= j;
796 }
797 break;
798 default:
799 if ( *term >= FUNCTION ) {
800 i = *term;
801 m = instruct;
802 while ( m < stopper ) {
803 if ( *m == -SNUMBER ) { NEXTARG(m) continue; }
804 if ( *m == -i && m[1] == -SNUMBER && (m+1) < stopper ) {
805 count += m[2]; m += 3;
806 }
807 NEXTARG(m)
808 }
809 if ( functions[i-FUNCTION].spec >= TENSORFUNCTION ) {
810 i = term[1] - FUNHEAD;
811 term += FUNHEAD;
812 while ( i > 0 ) {
813 if ( *term < 0 ) {
814 m = instruct;
815 while ( m < stopper ) {
816 if ( *m == -SNUMBER ) { NEXTARG(m) continue; }
817 if ( ( *m == -VECTOR || *m == -INDEX
818 || *m == -MINVECTOR ) && m[1] == *term &&
819 m[2] == -SNUMBER && (m+2) < stopper ) {
820 count += m[3]; m += 4;
821 }
822 else { NEXTARG(m) }
823 }
824 }
825 term++;
826 i--;
827 }
828 }
829 else {
830 r = term + term[1];
831 term += FUNHEAD;
832 while ( term < r ) {
833 if ( ( *term == -INDEX || *term == -VECTOR
834 || *term == -MINVECTOR ) && term[1] < MINSPEC ) {
835 m = instruct;
836 while ( m < stopper ) {
837 if ( *m == -SNUMBER ) { NEXTARG(m) continue; }
838 if ( *m == -VECTOR && m[1] == term[1]
839 && m[2] == -SNUMBER && (m+2) < stopper ) {
840 count += m[3];
841 m += 4;
842 }
843 else { NEXTARG(m) }
844 }
845 term += 2;
846 }
847 else { NEXTARG(term) }
848 }
849 }
850 break;
851 }
852 else {
853 term += term[1];
854 }
855 break;
856 }
857 }
858 return(count);
859}
860
861/*
862 #] CountFun :
863 #] Count :
864 #[ DimensionSubterm :
865*/
866
867WORD DimensionSubterm(WORD *subterm)
868{
869 WORD *r, *rstop, dim, i;
870 LONG x = 0;
871 rstop = subterm + subterm[1];
872 if ( *subterm == SYMBOL ) {
873 r = subterm + 2;
874 while ( r < rstop ) {
875 if ( *r <= NumSymbols && *r > -MAXPOWER ) {
876 dim = symbols[*r].dimension;
877 if ( dim == MAXPOSITIVE ) goto undefined;
878 x += dim * r[1];
879 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
880 r += 2;
881 }
882 else if ( *r <= MAXVARIABLES ) {
883/*
884 Here we have an extra symbol. Store dimension in the compiler buffer
885*/
886 i = MAXVARIABLES - *r;
887 dim = cbuf[AM.sbufnum].dimension[i];
888 if ( dim == MAXPOSITIVE ) goto undefined;
889 if ( dim == -MAXPOSITIVE ) goto outofrange;
890 x += dim * r[1];
891 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
892 r += 2;
893 }
894 else { r += 2; }
895 }
896 }
897 else if ( *subterm == DOTPRODUCT ) {
898 r = subterm + 2;
899 while ( r < rstop ) {
900 dim = vectors[*r-AM.OffsetVector].dimension;
901 if ( dim == MAXPOSITIVE ) goto undefined;
902 x += dim * r[2];
903 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
904 dim = vectors[r[1]-AM.OffsetVector].dimension;
905 if ( dim == MAXPOSITIVE ) goto undefined;
906 x += dim * r[2];
907 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
908 r += 3;
909 }
910 }
911 else if ( *subterm == VECTOR ) {
912 r = subterm + 2;
913 while ( r < rstop ) {
914 dim = vectors[*r-AM.OffsetVector].dimension;
915 if ( dim == MAXPOSITIVE ) goto undefined;
916 x += dim;
917 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
918 r += 2;
919 }
920 }
921 else if ( *subterm == INDEX ) {
922 r = subterm + 2;
923 while ( r < rstop ) {
924 if ( *r < 0 ) {
925 dim = vectors[*r-AM.OffsetVector].dimension;
926 if ( dim == MAXPOSITIVE ) goto undefined;
927 x += dim;
928 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
929 }
930 r++;
931 }
932 }
933 else if ( *subterm >= FUNCTION ) {
934 dim = functions[*subterm-FUNCTION].dimension;
935 if ( dim == MAXPOSITIVE ) goto undefined;
936 x += dim;
937 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
938 if ( functions[*subterm-FUNCTION].spec > 0 ) { /* tensor */
939 r = subterm + FUNHEAD;
940 while ( r < rstop ) {
941 if ( *r < MINSPEC ) {
942 dim = vectors[*r-AM.OffsetVector].dimension;
943 if ( dim == MAXPOSITIVE ) goto undefined;
944 x += dim;
945 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
946 }
947 r++;
948 }
949 }
950 }
951 return((WORD)x);
952undefined:
953 return((WORD)MAXPOSITIVE);
954outofrange:
955 return(-(WORD)MAXPOSITIVE);
956}
957
958/*
959 #] DimensionSubterm :
960 #[ DimensionTerm :
961
962 Returns the dimension of the given term.
963 If there is any variable of which the dimension is not defined
964 we return the code for undefined which is MAXPOSITIVE
965 When the value is out of range we return -MAXPOSITIVE
966*/
967
968WORD DimensionTerm(WORD *term)
969{
970 WORD *t, *tstop, dim;
971 LONG x = 0;
972 tstop = term + *term; tstop -= ABS(tstop[-1]);
973 t = term+1;
974 while ( t < tstop ) {
975 dim = DimensionSubterm(t);
976 if ( dim == MAXPOSITIVE ) goto undefined;
977 if ( dim == -MAXPOSITIVE ) goto outofrange;
978 x += dim;
979 if ( x >= MAXPOSITIVE || x <= -MAXPOSITIVE ) goto outofrange;
980 t += t[1];
981 }
982 return((WORD)x);
983undefined:
984 return((WORD)MAXPOSITIVE);
985outofrange:
986 return(-(WORD)MAXPOSITIVE);
987}
988
989/*
990 #] DimensionTerm :
991 #[ DimensionExpression :
992
993 Returns the dimension of the given expression.
994 If there is any variable of which the dimension is not defined
995 we return the code for undefined which is MAXPOSITIVE
996 When the value is out of range we return -MAXPOSITIVE
997 When the value is not consistent we return -MAXPOSITIVE.
998*/
999
1000WORD DimensionExpression(PHEAD WORD *expr)
1001{
1002 WORD dim, *term, *old, x = 0;
1003 int first = 1;
1004 term = expr;
1005 while ( *term ) {
1006 dim = DimensionTerm(term);
1007 if ( dim == MAXPOSITIVE ) goto undefined;
1008 if ( dim == -MAXPOSITIVE ) goto outofrange;
1009 if ( first ) { x = dim; }
1010 else if ( x != dim ) {
1011 old = AN.currentTerm;
1012 MLOCK(ErrorMessageLock);
1013 MesPrint("Dimension is not the same in the terms of the expression");
1014 term = expr;
1015 while ( *term ) {
1016 AN.currentTerm = term;
1017 MesPrint(" %T");
1018 }
1019 MUNLOCK(ErrorMessageLock);
1020 AN.currentTerm = old;
1021 return(-(WORD)MAXPOSITIVE);
1022 }
1023 term += *term;
1024 }
1025 return((WORD)x);
1026undefined:
1027 return((WORD)MAXPOSITIVE);
1028outofrange:
1029 old = AN.currentTerm;
1030 AN.currentTerm = term;
1031 MLOCK(ErrorMessageLock);
1032 MesPrint("Dimension out of range in %t in subexpression");
1033 MUNLOCK(ErrorMessageLock);
1034 AN.currentTerm = old;
1035 return(-(WORD)MAXPOSITIVE);
1036}
1037
1038/*
1039 #] DimensionExpression :
1040 #[ Multiply Term :
1041 #[ MultDo :
1042*/
1043
1044WORD MultDo(PHEAD WORD *term, WORD *pattern)
1045{
1046 GETBIDENTITY
1047 WORD *t, *r, i;
1048 t = term + *term;
1049 if ( pattern[2] > 0 ) { /* Left multiply */
1050 i = *term - 1;
1051 }
1052 else { /* Right multiply */
1053 i = ABS(t[-1]);
1054 }
1055 *term += SUBEXPSIZE;
1056 r = t + SUBEXPSIZE;
1057 do { *--r = *--t; } while ( --i > 0 );
1058 r = pattern + 3;
1059 i = r[1];
1060 while ( --i >= 0 ) *t++ = *r++;
1061 AT.WorkPointer = term + *term;
1062 return(0);
1063}
1064
1065/*
1066 #] MultDo :
1067 #] Multiply Term :
1068 #[ Try Term(s) :
1069 #[ TryDo :
1070*/
1071
1072WORD TryDo(PHEAD WORD *term, WORD *pattern, WORD level)
1073{
1074 GETBIDENTITY
1075 WORD *t, *r, *m, i, j;
1076 ReNumber(BHEAD term);
1077 Normalize(BHEAD term);
1078 m = r = term + *term;
1079 m++;
1080 i = pattern[2];
1081 t = pattern + 3;
1082 NCOPY(m,t,i)
1083 j = *term - 1;
1084 t = term + 1;
1085 NCOPY(m,t,j)
1086 *r = WORDDIF(m,r);
1087 AT.WorkPointer = m;
1088 if ( ( j = Normalize(BHEAD r) ) == 0 || j == 1 ) {
1089 if ( *r == 0 ) return(0);
1090 ReNumber(BHEAD r); Normalize(BHEAD r);
1091 if ( *r == 0 ) return(0);
1092 if ( ( i = CompareTerms(term,r,0) ) < 0 ) {
1093 *AN.RepPoint = 1;
1094 AR.expchanged = 1;
1095 return(Generator(BHEAD r,level));
1096 }
1097 if ( i == 0 && CompCoef(term,r) != 0 ) { return(0); }
1098 }
1099 AT.WorkPointer = r;
1100 return(Generator(BHEAD term,level));
1101}
1102
1103/*
1104 #] TryDo :
1105 #] Try Term(s) :
1106 #[ Distribute :
1107 #[ DoDistrib :
1108
1109 The routine that generates the terms ordered by a distrib_ function.
1110 The presence of a replaceable distrib_ function has been sensed
1111 in the routine TestSub and has been passed on to Generator.
1112 It is then Generator that calls this function in a way that is
1113 similar to calling the trace routines, except for that for the
1114 trace routines and the Levi-Civita tensors the arguments are put
1115 in temporary storage and here we leave them inside the term,
1116 because there is no knowing how long the field will be.
1117*/
1118
1119WORD DoDistrib(PHEAD WORD *term, WORD level)
1120{
1121 GETBIDENTITY
1122 WORD *t, *m, *r = 0, *stop, *tstop, *termout, *endhead, *starttail, *parms;
1123 WORD i, j, k, n, nn, ntype, fun1 = 0, fun2 = 0, typ1 = 0, typ2 = 0;
1124 WORD *arg, *oldwork, *mf, ktype = 0, atype = 0;
1125 WORD sgn, dirtyflag;
1126 AN.TeInFun = AR.TePos = 0;
1127 t = term;
1128 tstop = t + *t;
1129 stop = tstop - ABS(tstop[-1]);
1130 t++;
1131 while ( t < stop ) {
1132 r = t + t[1];
1133 if ( *t == DISTRIBUTION && t[FUNHEAD] == -SNUMBER
1134 && t[FUNHEAD+1] >= -2 && t[FUNHEAD+1] <= 2
1135 && t[FUNHEAD+2] == -SNUMBER
1136 && t[FUNHEAD+4] <= -FUNCTION
1137 && t[FUNHEAD+5] <= -FUNCTION ) {
1138 WORD *ttt = t+FUNHEAD+6, *tttstop = t+t[1];
1139 while ( ttt < tttstop ) {
1140 if ( *ttt == -DOLLAREXPRESSION ) break;
1141 NEXTARG(ttt);
1142 }
1143 if ( ttt >= tttstop ) {
1144 fun1 = -t[FUNHEAD+4];
1145 fun2 = -t[FUNHEAD+5];
1146 typ1 = functions[fun1-FUNCTION].spec;
1147 typ2 = functions[fun2-FUNCTION].spec;
1148 if ( typ1 > 0 || typ2 > 0 ) {
1149 m = t + FUNHEAD+6;
1150 r = t + t[1];
1151 while ( m < r ) {
1152 if ( *m != -INDEX && *m != -VECTOR && *m != -MINVECTOR )
1153 break;
1154 m += 2;
1155 }
1156 if ( m < r ) {
1157 MLOCK(ErrorMessageLock);
1158 MesPrint("Incompatible function types and arguments in distrib_");
1159 MUNLOCK(ErrorMessageLock);
1160 SETERROR(-1)
1161 }
1162 }
1163 break;
1164 }
1165 }
1166 t = r;
1167 }
1168 dirtyflag = t[2];
1169 ntype = t[FUNHEAD+1];
1170 n = t[FUNHEAD+3];
1171/*
1172 t points at the distrib_ function to be expanded.
1173 fun1,fun2 and typ1,typ2 are the two functions and their types.
1174 ntype indicates the action:
1175 0: Make all possible divisions: 2^nargs
1176 1: fun1 should get n arguments: nargs! / ( n! (nargs-n)! )
1177 2: fun2 should get n arguments: nargs! / ( n! (nargs-n)! )
1178 The distiction between 1 and two is for noncommuting objects.
1179 3: fun1 should get n arguments. Super symmetric option.
1180 4: fun2 idem
1181 The super symmetric option involves:
1182 a: arguments get sorted
1183 b: identical arguments are seen as such. Hence not all their
1184 distributions are taken into account. It is as if after the
1185 distrib there is a symmetrize fun1; symmetrize fun2;
1186 c: Hence if the occurrence of each argument is a,b,c,...
1187 and their occurrence in fun1 is a1,b1,c1,... and in fun2
1188 is a2,b2,c2,... then each term is generated (a1+a2)!/a1!/a2!
1189 (b1+b2)!/b1!/b2! (c1+c2)!/c1!/c2! ... times.
1190 d: We have to make an array of occurrences and positions.
1191 e: Then we sort the arguments indirectly.
1192 f: Next we generate the argument lists in the same way as we
1193 generate powers of expressions with binomials. Hence we need
1194 a third array to keep track of the `powers'
1195*/
1196 endhead = t;
1197 starttail = r;
1198 parms = m = t + FUNHEAD+6;
1199 i = 0;
1200 while ( m < r ) { /* Count arguments */
1201 i++;
1202 NEXTARG(m);
1203 }
1204 oldwork = AT.WorkPointer;
1205 arg = AT.WorkPointer + 1;
1206 arg[-1] = 0;
1207 termout = arg + i;
1208 switch ( ntype ) {
1209 case 0: ktype = 1; atype = n < 0 ? 1: 0; n = 0; break;
1210 case 1: ktype = 1; atype = 0; break;
1211 case 2: ktype = 0; atype = 0; break;
1212 case -1: ktype = 1; atype = 1; break;
1213 case -2: ktype = 0; atype = 1; break;
1214 }
1215 do {
1216/*
1217 All distributions with n elements. We generate the array arg with
1218 all possible 1 and 0 patterns. 1 means in fun1 and 0 means in fun2.
1219*/
1220 if ( n > i ) return(0); /* 0 elements */
1221
1222 for ( j = 0; j < n; j++ ) arg[j] = 1;
1223 for ( j = n; j < i; j++ ) arg[j] = 0;
1224 for(;;) {
1225 sgn = 0;
1226 t = term;
1227 m = termout;
1228 while ( t < endhead ) *m++ = *t++;
1229 mf = m;
1230 *m++ = fun1;
1231 *m++ = FUNHEAD;
1232 *m++ = dirtyflag;
1233#if FUNHEAD > 3
1234 k = FUNHEAD -3;
1235 while ( k-- > 0 ) *m++ = 0;
1236#endif
1237 r = parms;
1238 for ( k = 0; k < i; k++ ) {
1239 if ( arg[k] == ktype ) {
1240 if ( *r <= -FUNCTION ) *m++ = *r++;
1241 else if ( *r < 0 ) {
1242 if ( typ1 > 0 ) {
1243 if ( *r == -MINVECTOR ) sgn ^= 1;
1244 r++;
1245 *m++ = *r++;
1246 }
1247 else { *m++ = *r++; *m++ = *r++; }
1248 }
1249 else {
1250 nn = *r;
1251 NCOPY(m,r,nn);
1252 }
1253 }
1254 else { NEXTARG(r) }
1255 }
1256 mf[1] = WORDDIF(m,mf);
1257 mf = m;
1258 *m++ = fun2;
1259 *m++ = FUNHEAD;
1260 *m++ = dirtyflag;
1261#if FUNHEAD > 3
1262 k = FUNHEAD -3;
1263 while ( k-- > 0 ) *m++ = 0;
1264#endif
1265 r = parms;
1266 for ( k = 0; k < i; k++ ) {
1267 if ( arg[k] != ktype ) {
1268 if ( *r <= -FUNCTION ) *m++ = *r++;
1269 else if ( *r < 0 ) {
1270 if ( typ2 > 0 ) {
1271 if ( *r == -MINVECTOR ) sgn ^= 1;
1272 r++;
1273 *m++ = *r++;
1274 }
1275 else { *m++ = *r++; *m++ = *r++; }
1276 }
1277 else {
1278 nn = *r;
1279 NCOPY(m,r,nn);
1280 }
1281 }
1282 else { NEXTARG(r) }
1283 }
1284 mf[1] = WORDDIF(m,mf);
1285#ifndef NUOVO
1286 if ( atype == 0 ) {
1287 WORD k1,k2;
1288 for ( k = 0; k < i-1; k++ ) {
1289 if ( arg[k] == 0 ) continue;
1290 k1 = 1; k2 = k;
1291 while ( k < i-1 && EqualArg(parms,k,k+1) ) { k++; k1++; }
1292 while ( k2 <= k && arg[k2] == 1 ) k2++;
1293 k2 = k-k2+1;
1294/*
1295 Now we need k1!/(k2! (k1-k2)!)
1296*/
1297 if ( k2 != k1 && k2 != 0 ) {
1298 if ( GetBinom((UWORD *)m+3,m+2,k1,k2) ) {
1299 MLOCK(ErrorMessageLock);
1300 MesCall("DoDistrib");
1301 MUNLOCK(ErrorMessageLock);
1302 SETERROR(-1)
1303 }
1304 m[1] = ( m[2] < 0 ? -m[2]: m[2] ) + 3;
1305 *m = LNUMBER;
1306 m += m[1];
1307 }
1308 }
1309 }
1310#endif
1311 r = starttail;
1312 while ( r < tstop ) *m++ = *r++;
1313
1314 if ( atype ) { /* antisymmetric field */
1315 k = n;
1316 nn = 0;
1317 for ( j = 0; j < i && k > 0; j++ ) {
1318 if ( arg[j] == 1 ) k--;
1319 else nn += k;
1320 }
1321 sgn ^= nn & 1;
1322 }
1323
1324 if ( sgn ) m[-1] = -m[-1];
1325 *termout = WORDDIF(m,termout);
1326 AT.WorkPointer = m;
1327 if ( AT.WorkPointer > AT.WorkTop ) {
1328 MLOCK(ErrorMessageLock);
1329 MesWork();
1330 MUNLOCK(ErrorMessageLock);
1331 return(-1);
1332 }
1333 *AN.RepPoint = 1;
1334 AR.expchanged = 1;
1335 if ( Generator(BHEAD termout,level) ) Terminate(-1);
1336#ifndef NUOVO
1337 {
1338 WORD k1;
1339 j = i - 1;
1340 k = 0;
1341redok: while ( arg[j] == 1 && j >= 0 ) { j--; k++; }
1342 while ( arg[j] == 0 && j >= 0 ) j--;
1343 if ( j < 0 ) break;
1344 k1 = j;
1345 arg[j] = 0;
1346 while ( !atype && EqualArg(parms,j,j+1) ) {
1347 j++;
1348 if ( j >= i - k - 1 ) { j = k1; k++; goto redok; }
1349 arg[j] = 0;
1350 }
1351 while ( k >= 0 ) { j++; arg[j] = 1; k--; }
1352 j++;
1353 while ( j < i ) { arg[j] = 0; j++; }
1354 }
1355#else
1356 j = i - 1;
1357 k = 0;
1358 while ( arg[j] == 1 && j >= 0 ) { j--; k++; }
1359 while ( arg[j] == 0 && j >= 0 ) j--;
1360 if ( j < 0 ) break;
1361 arg[j] = 0;
1362 while ( k >= 0 ) { j++; arg[j] = 1; k--; }
1363 j++;
1364 while ( j < i ) { arg[j] = 0; j++; }
1365#endif
1366 }
1367 } while ( ntype == 0 && ++n <= i );
1368 AT.WorkPointer = oldwork;
1369 return(0);
1370}
1371
1372/*
1373 #] DoDistrib :
1374 #[ EqualArg :
1375
1376 Returns 1 if the arguments in the field are identical.
1377*/
1378
1379WORD EqualArg(WORD *parms, WORD num1, WORD num2)
1380{
1381 WORD *t1, *t2;
1382 WORD i;
1383 t1 = parms;
1384 while ( --num1 >= 0 ) { NEXTARG(t1); }
1385 t2 = parms;
1386 while ( --num2 >= 0 ) { NEXTARG(t2); }
1387 if ( *t1 != *t2 ) return(0);
1388 if ( *t1 < 0 ) {
1389 if ( *t1 <= -FUNCTION || t1[1] == t2[1] ) return(1);
1390 return(0);
1391 }
1392 i = *t1;
1393 while ( --i >= 0 ) {
1394 if ( *t1 != *t2 ) return(0);
1395 t1++; t2++;
1396 }
1397 return(1);
1398}
1399
1400/*
1401 #] EqualArg :
1402 #[ DoDelta3 :
1403*/
1404
1405WORD DoDelta3(PHEAD WORD *term, WORD level)
1406{
1407 GETBIDENTITY
1408 WORD *t, *m, *m1, *m2, *stopper, *tstop, *termout, *dels, *taken;
1409 WORD *ic, *jc, *factors;
1410 WORD num, num2, i, j, k, knum, a;
1411 AN.TeInFun = AR.TePos = 0;
1412 tstop = term + *term;
1413 stopper = tstop - ABS(tstop[-1]);
1414 t = term+1;
1415 while ( ( *t != DELTA3 || ((t[1]-FUNHEAD) & 1 ) != 0 ) && t < stopper )
1416 t += t[1];
1417 if ( t >= stopper ) {
1418 MLOCK(ErrorMessageLock);
1419 MesPrint("Internal error with dd_ function");
1420 MUNLOCK(ErrorMessageLock);
1421 Terminate(-1);
1422 }
1423 m1 = t; m2 = t + t[1];
1424 num = t[1] - FUNHEAD;
1425 if ( num == 0 ) {
1426 termout = t = AT.WorkPointer;
1427 m = term;
1428 while ( m < m1 ) *t++ = *m++;
1429 m = m2; while ( m < tstop ) *t++ = *m++;
1430 *termout = WORDDIF(t,termout);
1431 AT.WorkPointer = t;
1432 *AN.RepPoint = 1;
1433 AR.expchanged = 1;
1434 if ( Generator(BHEAD termout,level) ) {
1435 MLOCK(ErrorMessageLock);
1436 MesCall("Do dd_");
1437 MUNLOCK(ErrorMessageLock);
1438 SETERROR(-1)
1439 }
1440 AT.WorkPointer = termout;
1441 return(0);
1442 }
1443 t += FUNHEAD;
1444/*
1445 Step 1: sort the arguments
1446*/
1447 for ( i = 1; i < num; i++ ) {
1448 if ( t[i] < t[i-1] ) {
1449 a = t[i]; t[i] = t[i-1]; t[i-1] = a;
1450 j = i - 1;
1451 while ( j > 0 ) {
1452 if ( t[j] >= t[j-1] ) break;
1453 a = t[j]; t[j] = t[j-1]; t[j-1] = a;
1454 j--;
1455 }
1456 }
1457 }
1458/*
1459 Step 2: Order them by occurrence
1460 In 'taken' we have the array with the number of occurrences.
1461 in 'dels' is the type of object.
1462*/
1463 m = taken = AT.WorkPointer;
1464 for ( i = 0; i < num; i++ ) *m++ = 0;
1465 dels = m; knum = 0;
1466 for ( i = 0; i < num; knum++ ) {
1467 *m++ = t[i]; i++; taken[knum] = 1;
1468 while ( i < num ) {
1469 if ( t[i] != t[i-1] ) break;
1470 i++; (taken[knum])++;
1471 }
1472 }
1473 for ( i = 0; i < knum; i++ ) *m++ = taken[i];
1474 ic = m; num2 = num/2;
1475 jc = ic + num2;
1476 factors = jc + num2;
1477 termout = factors + num2;
1478/*
1479 The recursion has num/2 steps
1480*/
1481 k = 0;
1482 while ( k >= 0 ) {
1483 if ( k >= num2 ) {
1484 t = termout; m = term;
1485 while ( m < m1 ) *t++ = *m++;
1486 *t++ = DELTA; *t++ = num+2;
1487 for ( i = 0; i < num2; i++ ) {
1488 *t++ = dels[ic[i]]; *t++ = dels[jc[i]];
1489 }
1490 for ( i = 0; i < num2; i++ ) {
1491 if ( ic[i] == jc[i] ) {
1492 j = 1;
1493 while ( i < num2-1 && ic[i] == ic[i+1] && ic[i] == jc[i+1] )
1494 { i++; j++; }
1495 for ( a = 1; a < j; a++ ) {
1496 *t++ = SNUMBER; *t++ = 4; *t++ = 2*a+1; *t++ = 1;
1497 }
1498 for ( a = 0; a+1+i < num2; a++ ) {
1499 if ( ic[a+i] != ic[a+i+1] ) break;
1500 }
1501 if ( a > 0 ) {
1502 if ( GetBinom((UWORD *)(t+3),t+2,2*j+a,a) ) {
1503 MLOCK(ErrorMessageLock);
1504 MesCall("Do dd_");
1505 MUNLOCK(ErrorMessageLock);
1506 SETERROR(-1)
1507 }
1508 t[1] = ( t[2] < 0 ? -t[2]: t[2] ) + 3;
1509 *t = LNUMBER;
1510 t += t[1];
1511 }
1512 }
1513 else if ( factors[i] != 1 ) {
1514 *t++ = SNUMBER; *t++ = 4; *t++ = factors[i]; *t++ = 1;
1515 }
1516 }
1517 for ( i = 0; i < num2-1; i++ ) {
1518 if ( ic[i] == jc[i] ) continue;
1519 j = 1;
1520 while ( i < num2-1 && jc[i] == jc[i+1] && ic[i] == ic[i+1] ) {
1521 i++; j++;
1522 }
1523 for ( a = 0; a+i < num2-1; a++ ) {
1524 if ( ic[i+a] != ic[i+a+1] ) break;
1525 }
1526 if ( a > 0 ) {
1527 if ( GetBinom((UWORD *)(t+3),t+2,j+a,a) ) {
1528 MLOCK(ErrorMessageLock);
1529 MesCall("Do dd_");
1530 MUNLOCK(ErrorMessageLock);
1531 SETERROR(-1)
1532 }
1533 t[1] = ( t[2] < 0 ? -t[2]: t[2] ) + 3;
1534 *t = LNUMBER;
1535 t += t[1];
1536 }
1537 }
1538 m = m2;
1539 while ( m < tstop ) *t++ = *m++;
1540 *termout = WORDDIF(t,termout);
1541 AT.WorkPointer = t;
1542 *AN.RepPoint = 1;
1543 AR.expchanged = 1;
1544 if ( Generator(BHEAD termout,level) ) {
1545 MLOCK(ErrorMessageLock);
1546 MesCall("Do dd_");
1547 MUNLOCK(ErrorMessageLock);
1548 SETERROR(-1)
1549 }
1550 k--;
1551 if ( k >= 0 ) goto nextj;
1552 else break;
1553 }
1554 for ( ic[k] = 0; ic[k] < knum; ic[k]++ ) {
1555 if ( taken[ic[k]] > 0 ) break;
1556 }
1557 if ( k > 0 && ic[k-1] == ic[k] ) jc[k] = jc[k-1];
1558 else jc[k] = ic[k];
1559 for ( ; jc[k] < knum; jc[k]++ ) {
1560 if ( taken[jc[k]] <= 0 ) continue;
1561 if ( ic[k] == jc[k] ) {
1562 if ( taken[jc[k]] <= 1 ) continue;
1563/*
1564 factors[k] = taken[ic[k]];
1565 if ( ( factors[k] & 1 ) == 0 ) (factors[k])--;
1566*/
1567 taken[ic[k]] -= 2;
1568 }
1569 else {
1570 factors[k] = taken[jc[k]];
1571 (taken[ic[k]])--; (taken[jc[k]])--;
1572 }
1573 k++;
1574 goto nextk; /* This is the simulated recursion */
1575nextj:;
1576 (taken[ic[k]])++; (taken[jc[k]])++;
1577 }
1578 k--;
1579 if ( k >= 0 ) goto nextj;
1580nextk:;
1581 }
1582 AT.WorkPointer = taken;
1583 return(0);
1584}
1585
1586/*
1587 #] DoDelta3 :
1588 #[ TestPartitions :
1589
1590 Checks whether the function in tfun is a partitions_ function
1591 that can be expanded. If it can a number of relevant objects is
1592 inside the struct parti.
1593 This test is not entirely trivial because there are many restrictions
1594 w.r.t. the arguments.
1595 Syntax (still to be implemented)
1596 partitions_(number_of_partition_entries,[function,number,]^nope,arguments)
1597 [function,number,]: can be
1598 f,3 for a partition of 3 arguments
1599 f,0 for the remaining arguments (should be last)
1600 num1,f,num2 with num1 effectively a number of partitions but this
1601 counts as num1 entries.
1602 0,f,num2: all partitions have num2 arguments. No number of partition
1603 entries needed. If num2 does not divide the number of
1604 arguments there will be no action.
1605*/
1606
1607WORD TestPartitions(WORD *tfun, PARTI *parti)
1608{
1609 WORD *tnext = tfun + tfun[1];
1610 WORD *t, *tt;
1611 WORD argcount = 0, sum = 0, i, ipart, argremain;
1612 WORD tensorflag = 0;
1613 parti->psize = parti->nfun = parti->args = parti->nargs = 0;
1614 parti->numargs = parti->numpart = parti->where = 0;
1615 tt = t = tfun + FUNHEAD;
1616 while ( t < tnext ) { argcount++; NEXTARG(t); }
1617 if ( argcount < 1 ) goto No;
1618 t = tt;
1619 if ( *t != -SNUMBER ) goto No;
1620 if ( t[1] == 0 ) {
1621 t += 2;
1622 if ( *t <= -FUNCTION && t[1] == -SNUMBER && t[2] > 0 ) {
1623 if ( functions[-*t-FUNCTION].spec > 0 ) tensorflag = 1;
1624 if ( argcount-3 < 0 ) goto No;
1625 if ( ( (argcount-3) % t[2] ) != 0 ) goto No;
1626 }
1627 else goto No;
1628 parti->numpart = (argcount-3)/t[2];
1629 parti->numargs = argcount - 3;
1630 parti->psize = (WORD *)Malloc1((parti->numpart*2+parti->numargs*2+2)
1631 *sizeof(WORD),"partitions");
1632 parti->nfun = parti->psize + parti->numpart;
1633 parti->args = parti->nfun + parti->numpart;
1634 parti->nargs = parti->args + parti->numargs;
1635 for ( i = 0; i < parti->numpart; i++ ) {
1636 parti->psize[i] = t[2];
1637 parti->nfun[i] = -t[0];
1638 }
1639 t += 3;
1640 }
1641 else if ( t[1] > 0 ) { /* Number of partitions */
1642/*
1643 We can have sequences of function,number for one partition
1644 or number1,function,number2 for number1 partitions of size number2.
1645 The last partition can have number=0. It must be a single partition
1646 and it will take all remaining arguments.
1647 If any of the functions is a tensor, all arguments must be either
1648 vector or index.
1649*/
1650 parti->numpart = t[1]; t += 2;
1651 ipart = sum = 0; argremain = argcount - 1;
1652/*
1653 At this point is seems better to make an allocation already that
1654 may be too big. The alternative is having to pass this code twice.
1655*/
1656 parti->psize = (WORD *)Malloc1((argcount*4+2)*sizeof(WORD),"partitions");
1657 parti->nfun = parti->psize+argcount;
1658 parti->args = parti->nfun+argcount;
1659 parti->nargs = parti->args+argcount;
1660 while ( ipart < parti->numpart ) {
1661 if ( *t <= -FUNCTION && t[1] == -SNUMBER && t[2] >= 0 ) {
1662 if ( functions[-*t-FUNCTION].spec > 0 ) tensorflag = 1;
1663 if ( t[2] == 0 ) {
1664 if ( ipart+1 != parti->numpart ) goto WhatAPity;
1665 argremain -= 2;
1666 parti->nfun[ipart] = -*t;
1667 parti->psize[ipart++] = argremain-sum;
1668 ipart++;
1669 sum = argremain;
1670 }
1671 else {
1672 parti->nfun[ipart] = -*t;
1673 parti->psize[ipart++] = t[2];
1674 argremain -= 2;
1675 sum += t[2];
1676 }
1677 t += 3;
1678 }
1679 else if ( *t == -SNUMBER && t[1] > 0 && ipart+t[1] <= parti->numpart
1680 && t[2] <= -FUNCTION && t[3] == -SNUMBER && t[4] > 0 ) {
1681 if ( functions[-t[2]-FUNCTION].spec > 0 ) tensorflag = 1;
1682 argremain -= 3;
1683 for ( i = 0; i < t[1]; i++ ) {
1684 parti->nfun[ipart] = -t[2];
1685 parti->psize[ipart++] = t[4];
1686 sum += t[4];
1687 }
1688 if ( sum > argremain ) goto WhatAPity;
1689 t += 5;
1690 }
1691 else goto WhatAPity;
1692 }
1693 if ( sum != argremain ) goto WhatAPity;
1694 parti->numargs = argremain;
1695 }
1696 else goto No;
1697/*
1698 Now load the offsets of the arguments and check if needed whether OK with tensor
1699*/
1700 for ( i = 0; i < parti->numargs; i++ ) {
1701 parti->args[i] = t - tfun;
1702 if ( tensorflag && ( *t != -VECTOR && *t != -INDEX ) ) goto WhatAPity;
1703 NEXTARG(t);
1704 }
1705 return(1);
1706WhatAPity:
1707 M_free(parti->psize,"partitions");
1708 parti->psize = parti->nfun = parti->args = parti->nargs = 0;
1709 parti->numargs = parti->numpart = parti->where = 0;
1710No:
1711 return(0);
1712}
1713
1714/*
1715 #] TestPartitions :
1716 #[ DoPartitions :
1717
1718 As we have only one AT.partitions we need to copy it locally
1719 if we keep needing it.
1720*/
1721
1722WORD DoPartitions(PHEAD WORD *term, WORD level)
1723{
1724 WORD x, i, j, im, *fun, ndiff, siz, tensorflag = 0;
1725 PARTI part = AT.partitions;
1726 WORD *array, **j3, **j3fill, **j3where;
1727 WORD a, pfill, *j2, *j2fill, j3size, ncoeff, ncoeffnum, nfac, ncoeff2, ncoeff3, n;
1728 UWORD *coeff, *coeffnum, *cfac, *coeff2, *coeff3, *c;
1729 /* Make AT.partitions ready for future use (if there is another function) */
1730 AT.partitions.psize = AT.partitions.nfun = AT.partitions.args = AT.partitions.nargs = 0;
1731 AT.partitions.numargs = AT.partitions.numpart = AT.partitions.where = 0;
1732/*
1733 Start with bubble sorting the list of arguments. And the list of partitions.
1734*/
1735 fun = term + part.where;
1736 if ( functions[*fun-FUNCTION].spec ) tensorflag = 1;
1737 for ( i = 1; i < part.numargs; i++ ) {
1738 for ( j = i-1; j >= 0; j-- ) {
1739 if ( CompArg(fun+part.args[j+1],fun+part.args[j]) >= 0 ) break;
1740 x = part.args[j+1]; part.args[j+1] = part.args[j]; part.args[j] = x;
1741 }
1742 }
1743 for ( i = 1; i < part.numpart; i++ ) {
1744 for ( j = i-1; j >= 0; j-- ) {
1745 if ( part.psize[j+1] < part.psize[j] ) break;
1746 if ( part.psize[j+1] == part.psize[j] && part.nfun[j+1] <= part.nfun[j] ) break;
1747 x = part.psize[j+1]; part.psize[j+1] = part.psize[j]; part.psize[j] = x;
1748 x = part.nfun[j+1]; part.nfun[j+1] = part.nfun[j]; part.nfun[j] = x;
1749 }
1750 }
1751/*
1752 Now we have the partitions sorted from high to low and the arguments
1753 have been sorted the regular way arguments are sorted in a symmetrize.
1754 The important thing is that identical arguments are adjacent.
1755 Assign the numbers (identical arguments have identical numbers).
1756*/
1757 ndiff = 1; part.nargs[0] = ndiff;
1758 for ( i = 1; i < part.numargs; i++ ) {
1759 if ( CompArg(fun+part.args[i],fun+part.args[i-1]) != 0 ) ndiff++;
1760 part.nargs[i] = ndiff;
1761 }
1762 part.nargs[part.numargs] = 0;
1763 coeffnum = NumberMalloc("partitionsn");
1764 coeff = NumberMalloc("partitions");
1765 coeff2 = NumberMalloc("partitions2");
1766 coeff3 = NumberMalloc("partitions3");
1767 cfac = NumberMalloc("partitions!");
1768 ncoeffnum = 1; coeffnum[0] = 1;
1769/*
1770 The numerator of the coefficient will be n1!*n2!*...*n(ndiff)!
1771 We compute it only once.
1772*/
1773 j = 0;
1774 for ( i = 1; i <= ndiff; i++ ) {
1775 n = 0;
1776 while ( part.nargs[j] == i ) { n++; j++; }
1777 if ( n > 1 ) { /* 1! needs no attention */
1778 if ( Factorial(BHEAD n, cfac, &nfac) ) Terminate(-1);
1779 if ( MulLong(coeffnum,ncoeffnum,cfac,nfac,coeff2,&ncoeff2) ) Terminate(-1);
1780 c = coeffnum; coeffnum = coeff2; coeff2 = c;
1781 n = ncoeffnum; ncoeffnum = ncoeff2; ncoeff2 = n;
1782 }
1783 }
1784/*
1785 Now comes the part where we have to make sure that
1786 a: we generate all partitions.
1787 b: we generate only different partitions.
1788 c: we get the proper combinatorics factor.
1789 Method:
1790 Suppose the largest partition needs n objects and there are m partitions.
1791 We allocate m arrays of n 'digits'. Make in the smaller partitions the
1792 appropriate leading digits zero.
1793 Divide the largest numbers (of the arguments) over the partitions as
1794 leftmost digits (after possible zeroes). The arrays, seen as numbers,
1795 should be such that each is less or equal to its left neighbour. Take the
1796 next largest numbers, etc. This generates unique partitions and all of
1797 them. Because we have a formula for the multiplicity, this should do it.
1798
1799 The general case. At a later stage we might put in a more economical
1800 version for special cases.
1801*/
1802 siz = part.psize[0];
1803 j3size = 2*(part.numpart+1)+2*(part.numargs+1);
1804 array = (WORD *)Malloc1((part.numpart+1)*siz*sizeof(WORD),"parts");
1805 j3 = (WORD **)Malloc1(j3size*sizeof(WORD *),"parts3");
1806 j2 = (WORD *)Malloc1((part.numpart+part.numargs+2)*sizeof(WORD),"parts2");
1807 j3fill = j3+(part.numpart+1);
1808 j3where = j3fill+(part.numpart+1);
1809 for ( i = 0; i < j3size; i++ ) j3[i] = 0;
1810 j2fill = j2+(part.numpart+1);
1811 for ( i = 0; i < part.numargs; i++ ) j2fill[i] = 0;
1812 for ( i = 0; i < part.numpart; i++ ) {
1813 j3[i] = array+i*siz;
1814 for ( j = 0; j < siz; j++ ) j3[i][j] = 0;
1815 j3fill[i] = j3[i]+(siz-part.psize[i]);
1816 j2[i] = part.psize[i]; /* Number of places still available */
1817 }
1818 j3[part.numpart] = array+part.numpart*siz;
1819 j2[part.numpart] = 0;
1820/*
1821 Now comes a complicated two-level recursion in a and pfill.
1822*/
1823 a = part.numargs-1;
1824 pfill = 0;
1825/*
1826 We start putting the last number in part.nargs in the first partition in array.
1827 For backtracking we need to know where we put this number. Hence j3where.
1828*/
1829 while ( a < part.numargs ) {
1830 while ( j2[pfill] <= 0 ) {
1831 pfill++;
1832 while ( pfill >= part.numpart ) { /* we have to pop */
1833 a++;
1834 if ( a >= part.numargs ) goto Done;
1835 pfill = j2fill[a];
1836 j2[pfill]++;
1837 j3where[a][0] = 0;
1838 j3fill[pfill]--;
1839 pfill++;
1840 }
1841 }
1842 j3where[a] = j3fill[pfill];
1843 *(j3fill[pfill])++ = part.nargs[a];
1844 j2[pfill]--; j2fill[a] = pfill;
1845/*
1846 Now test whether this is allowed.
1847*/
1848 if ( pfill > 0 && part.psize[pfill] == part.psize[pfill-1]
1849 && part.nfun[pfill] == part.nfun[pfill-1] ) { /* First check whether allowed */
1850 for ( im = 0; im < siz; im++ ) {
1851 if ( j3[pfill-1][im] < j3[pfill][im] ) break;
1852 if ( j3[pfill-1][im] > j3[pfill][im] ) im = siz;
1853 }
1854 if ( im < siz ) { /* not ordered. undo and raise pfill */
1855 pfill = j2fill[a];
1856 j2[pfill]++;
1857 j3where[a][0] = 0;
1858 j3fill[pfill]--;
1859 pfill++;
1860 continue; /* Note that j2[part.numpart] = 0 */
1861 }
1862 }
1863 a--;
1864 if ( a < 0 ) { /* Solution */
1865/*
1866 #[ Solution :
1867
1868 Now we compose the output term. The input term contains
1869 three parts: head, partitions_, tail.
1870 partitions_ starts at term+part.where.
1871 We first put the function parts and worry about the coefficient later.
1872*/
1873 WORD *t, *to, *twhere = term+part.where, *t2, *tend = term+*term, *termout;
1874 WORD num, jj, *targ, *tfun;
1875 t2 = twhere+twhere[1];
1876 to = termout = AT.WorkPointer;
1877 if ( termout + *term + part.numpart*FUNHEAD + AM.MaxTal >= AT.WorkTop ) {
1878 return(MesWork());
1879 }
1880 for ( i = 0; i < ncoeffnum; i++ ) coeff[i] = coeffnum[i];
1881 ncoeff = ncoeffnum;
1882 t = term; while ( t < twhere ) *to++ = *t++;
1883/*
1884 Now the partitions
1885*/
1886 for ( i = 0; i < part.numpart; i++ ) {
1887 tfun = to;
1888 *to++ = part.nfun[i]; to++; FILLFUN(to);
1889 for ( j = 1; j <= part.psize[i]; j++ ) {
1890 num = j3[i][siz-j]; /* now we need an argument with this number */
1891 for ( jj = num-1; jj < part.numargs; jj++ ) {
1892 if ( part.nargs[jj] == num ) break;
1893 }
1894 targ = part.args[jj]+twhere;
1895 if ( *targ < 0 ) {
1896 if ( tensorflag ) targ++;
1897 else if ( *targ > -FUNCTION ) *to++ = *targ++;
1898 *to++ = *targ++;
1899 }
1900 else { jj = *targ; NCOPY(to,targ,jj); }
1901 }
1902 tfun[1] = to - tfun;
1903 }
1904/*
1905 Now the denominators of the coefficient
1906 First identical functions/partitions
1907*/
1908 j = 1; n = 1;
1909 while ( j < part.numpart ) {
1910 for ( im = 0; im < siz; im++ ) {
1911 if ( part.nfun[j-1] != part.nfun[j] ) break;
1912 if ( j3[j-1][im] < j3[j][im] ) break;
1913 if ( j3[j-1][im] > j3[j][im] ) im = 2*siz+2;
1914 }
1915 if ( im == siz ) { n++; j++; continue; }
1916 if ( n > 1 ) {
1917div1: if ( Factorial(BHEAD n, cfac, &nfac) ) Terminate(-1);
1918 if ( DivLong(coeff,ncoeff,cfac,nfac,coeff2,&ncoeff2,coeff3,&ncoeff3) ) Terminate(-1);
1919 c = coeff; coeff = coeff2; coeff2 = c;
1920 n = ncoeff; ncoeff = ncoeff2; ncoeff2 = n;
1921 }
1922 n = 1; j++;
1923 }
1924 if ( n > 1 ) goto div1;
1925/*
1926 Now identical elements inside the partitions
1927*/
1928 for ( i = 0; i < part.numpart; i++ ) {
1929 j = 0; while ( j3[i][j] == 0 ) j++;
1930 n = 1; j++;
1931 while ( j < siz ) {
1932 if ( j3[i][j-1] == j3[i][j] ) { n++; j++; }
1933 else {
1934 if ( n > 1 ) {
1935div2: if ( Factorial(BHEAD n, cfac, &nfac) ) Terminate(-1);
1936 if ( DivLong(coeff,ncoeff,cfac,nfac,coeff2,&ncoeff2,coeff3,&ncoeff3) ) Terminate(-1);
1937 c = coeff; coeff = coeff2; coeff2 = c;
1938 n = ncoeff; ncoeff = ncoeff2; ncoeff2 = n;
1939 }
1940 n = 1; j++;
1941 }
1942 }
1943 if ( n > 1 ) goto div2;
1944 }
1945/*
1946 And put this inside the term. Normalize will take care of it.
1947*/
1948 if ( ncoeff != 1 || coeff[0] > 1 ) {
1949 if ( ncoeff == 1 && coeff[0] <= MAXPOSITIVE ) {
1950 *to++ = SNUMBER; *to++ = 4; *to++ = (WORD)(coeff[0]); *to++ = 1;
1951 }
1952 else {
1953 *to++ = LNUMBER; *to++ = ncoeff+3; *to++ = ncoeff;
1954 for ( i = 0; i < ncoeff; i++ ) *to++ = ((WORD *)coeff)[i];
1955 }
1956 }
1957/*
1958 And the tail
1959*/
1960 while ( t2 < tend ) *to++ = *t2++;
1961 *termout = to-termout;
1962 AT.WorkPointer = to;
1963 if ( Generator(BHEAD termout,level) ) Terminate(-1);
1964 AT.WorkPointer = termout;
1965/*
1966 #] Solution :
1967
1968 Now we can pop all a with the lowest value and one more.
1969*/
1970 a = 0;
1971 while ( part.nargs[a] == 1 ) {
1972 pfill = j2fill[a]; j2[pfill]++; j3where[a][0] = 0; j3fill[pfill]--; a++;
1973 }
1974 if ( a < part.numargs ) {
1975 pfill = j2fill[a]; j2[pfill]++; j3where[a][0] = 0; j3fill[pfill]--; a++;
1976 }
1977 a--;
1978 pfill++;
1979 }
1980 else if ( part.nargs[a] == part.nargs[a+1] ) {}
1981 else { pfill = 0; }
1982 }
1983Done:
1984 M_free(j2,"parts2");
1985 M_free(j3,"parts3");
1986 M_free(array,"parts");
1987 NumberFree(cfac,"partitions!");
1988 NumberFree(coeff3,"partitions3");
1989 NumberFree(coeff2,"partitions2");
1990 NumberFree(coeff,"partitions");
1991 NumberFree(coeffnum,"partitionsn");
1992 M_free(part.psize,"partitions");
1993 part.psize = part.nfun = part.args = part.nargs = 0;
1994 part.numargs = part.numpart = part.where = 0;
1995 return(0);
1996}
1997
1998/*
1999 #] DoPartitions :
2000 #] Distribute :
2001 #[ DoPermutations :
2002
2003 Routine replaces the function perm_(f,args) by occurrences of f with
2004 all permutations of the args. This should always fit!
2005*/
2006
2007WORD DoPermutations(PHEAD WORD *term, WORD level)
2008{
2009 PERMP perm;
2010 WORD *oldworkpointer = AT.WorkPointer, *termout = AT.WorkPointer;
2011 WORD *t, *tstop, *tt, *ttstop, odd = 0;
2012 WORD *args[MAXMATCH], nargs, i, first, skip, *to, *from;
2013/*
2014 Find function and count arguments. Check for odd/even
2015*/
2016 tstop = term+*term; tstop -= ABS(tstop[-1]);
2017 t = term+1;
2018 while ( t < tstop ) {
2019 if ( *t == PERMUTATIONS ) {
2020 if ( t[1] >= FUNHEAD+1 && t[FUNHEAD] <= -FUNCTION ) {
2021 odd = 0; skip = 1;
2022 }
2023 else if ( t[1] >= FUNHEAD+3 && t[FUNHEAD] == -SNUMBER && t[FUNHEAD+2] <= -FUNCTION ) {
2024 if ( t[FUNHEAD+1] % 2 == 1 ) odd = -1;
2025 else odd = 0;
2026 skip = 3;
2027 }
2028 else { t += t[1]; continue; }
2029 tt = t+FUNHEAD+skip; ttstop = t + t[1];
2030 nargs = 0;
2031 while ( tt < ttstop ) { NEXTARG(tt); nargs++; }
2032 tt = t+FUNHEAD+skip;
2033 if ( nargs > MAXMATCH ) {
2034 MLOCK(ErrorMessageLock);
2035 MesPrint("Too many arguments in function perm_. %d! is way too big",(WORD)MAXMATCH);
2036 MUNLOCK(ErrorMessageLock);
2037 SETERROR(-1)
2038 }
2039 i = 0;
2040 while ( tt < ttstop ) { args[i++] = tt; NEXTARG(tt); }
2041 perm.n = nargs;
2042 perm.sign = 0;
2043 perm.objects = args;
2044 first = 1;
2045 while ( (first = PermuteP(&perm,first) ) == 0 ) {
2046/*
2047 Compose the output term
2048*/
2049 to = termout; from = term;
2050 while ( from < t ) *to++ = *from++;
2051 *to++ = -t[FUNHEAD+skip-1];
2052 *to++ = t[1] - skip;
2053 for ( i = 2; i < FUNHEAD; i++ ) *to++ = t[i];
2054 for ( i = 0; i < nargs; i++ ) {
2055 from = args[i];
2056 COPY1ARG(to,from);
2057 }
2058 from = t+t[1];
2059 tstop = term + *term;
2060 while ( from < tstop ) *to++ = *from++;
2061 if ( odd && ( ( perm.sign & 1 ) != 0 ) ) to[-1] = -to[-1];
2062 *termout = to - termout;
2063 AT.WorkPointer = to;
2064 if ( Generator(BHEAD termout,level) ) Terminate(-1);
2065 AT.WorkPointer = oldworkpointer;
2066 }
2067 return(0);
2068 }
2069 t += t[1];
2070 }
2071 return(0);
2072}
2073
2074/*
2075 #] DoPermutations :
2076 #[ DoShuffle :
2077
2078 Merges the arguments of all occurrences of function fun into a
2079 single occurrence of fun. The opposite of Distrib_
2080 Syntax:
2081 Shuffle[,once|all],fun;
2082 Shuffle[,once|all],$fun;
2083 The expansion of the dollar should give a single function.
2084 The dollar is indicated as usual with a negative value.
2085 option = 1 (once): generate identical results only once
2086 option = 0 (all): generate identical results with combinatorics (default)
2087*/
2088
2089/*
2090 We use the Shuffle routine which has a large amount of combinatorics.
2091 It doesn't have grouped combinatorics as in (0,1,2)*(0,1,3) where the
2092 groups (0,1) also cause double terms.
2093*/
2094
2095WORD DoShuffle(WORD *term, WORD level, WORD fun, WORD option)
2096{
2097 GETIDENTITY
2098 SHvariables SHback, *SH = &(AN.SHvar);
2099 WORD *t1, *t2, *tstop, ncoef, n = fun, *to, *from;
2100 int i, error;
2101 LONG k;
2102 UWORD *newcombi;
2103
2104 if ( n < 0 ) {
2105 if ( ( n = DolToFunction(BHEAD -n) ) == 0 ) {
2106 MLOCK(ErrorMessageLock);
2107 MesPrint("$-variable in merge statement did not evaluate to a function.");
2108 MUNLOCK(ErrorMessageLock);
2109 return(1);
2110 }
2111 }
2112 if ( AT.WorkPointer + 3*(*term) + AM.MaxTal > AT.WorkTop ) {
2113 MLOCK(ErrorMessageLock);
2114 MesWork();
2115 MUNLOCK(ErrorMessageLock);
2116 return(-1);
2117 }
2118
2119 tstop = term + *term;
2120 ncoef = tstop[-1];
2121 tstop -= ABS(ncoef);
2122 t1 = term + 1;
2123 while ( t1 < tstop ) {
2124 if ( ( *t1 == n ) && ( t1+t1[1] < tstop ) && ( t1[1] > FUNHEAD ) ) {
2125 t2 = t1 + t1[1];
2126 if ( t2 >= tstop ) {
2127 return(Generator(BHEAD term,level));
2128 }
2129 while ( t2 < tstop ) {
2130 if ( ( *t2 == n ) && ( t2[1] > FUNHEAD ) ) break;
2131 t2 += t2[1];
2132 }
2133 if ( t2 < tstop ) break;
2134 }
2135 t1 += t1[1];
2136 }
2137 if ( t1 >= tstop ) {
2138 return(Generator(BHEAD term,level));
2139 }
2140 *AN.RepPoint = 1;
2141/*
2142 Now we have two occurrences of the function.
2143 Back up all relevant variables and load all the stuff that needs to be
2144 passed on.
2145*/
2146 SHback = AN.SHvar;
2147 SH->finishuf = &FinishShuffle;
2148 SH->do_uffle = &DoShuffle;
2149 SH->outterm = AT.WorkPointer;
2150 AT.WorkPointer += *term;
2151 SH->stop1 = t1 + t1[1];
2152 SH->stop2 = t2 + t2[1];
2153 SH->thefunction = n;
2154 SH->option = option;
2155 SH->level = level;
2156 SH->incoef = tstop;
2157 SH->nincoef = ncoef;
2158
2159 if ( AN.SHcombi == 0 || AN.SHcombisize == 0 ) {
2160 AN.SHcombisize = 200;
2161 AN.SHcombi = (UWORD *)Malloc1(AN.SHcombisize*sizeof(UWORD),"AN.SHcombi");
2162 SH->combilast = 0;
2163 SHback.combilast = 0;
2164 }
2165 else {
2166 SH->combilast += AN.SHcombi[SH->combilast]+1;
2167 if ( SH->combilast >= AN.SHcombisize - 100 ) {
2168 newcombi = (UWORD *)Malloc1(2*AN.SHcombisize*sizeof(UWORD),"AN.SHcombi");
2169 for ( k = 0; k < AN.SHcombisize; k++ ) newcombi[k] = AN.SHcombi[k];
2170 M_free(AN.SHcombi,"AN.SHcombi");
2171 AN.SHcombi = newcombi;
2172 AN.SHcombisize *= 2;
2173 }
2174 }
2175 AN.SHcombi[SH->combilast] = 1;
2176 AN.SHcombi[SH->combilast+1] = 1;
2177
2178 i = t1-term; to = SH->outterm; from = term;
2179 NCOPY(to,from,i)
2180 SH->outfun = to;
2181 for ( i = 0; i < FUNHEAD; i++ ) { *to++ = t1[i]; }
2182
2183 error = Shuffle(t1+FUNHEAD,t2+FUNHEAD,to);
2184
2185 AT.WorkPointer = SH->outterm;
2186 AN.SHvar = SHback;
2187 if ( error ) {
2188 MesCall("DoShuffle");
2189 return(-1);
2190 }
2191 return(0);
2192}
2193
2194/*
2195 #] DoShuffle :
2196 #[ Shuffle :
2197
2198 How to make shuffles:
2199
2200 We have two lists of arguments. We have to make a single
2201 shuffle of them. All combinations. Doubles should have as
2202 much as possible a combinatorics factor. Sometimes this is
2203 very difficult as in:
2204 (0,1,2)x(0,1,3) = -> (0,1) is a repeated pattern and the
2205 factor on that is difficult
2206 Simple way: (without combinatorics)
2207 repeat id f0(?c)*f(x1?,?a)*f(x2?,?b) =
2208 +f0(?c,x1)*f(?a)*f(x2,?b)
2209 +f0(?c,x2)*f(x1,?a)*f(?b);
2210 Refinement:
2211 if ( x1 == x2 ) check how many more there are of the same.
2212 --> (n1,x) and (n2,x)
2213 id f0(?c)*f1((n1,x),?b)*f2((n2,x),?c) =
2214 +binom_(n1+n2,n1)*f0(?c,(n1+n2,x))*f1(?a)*f2(?b)
2215 +sum_(j,0,n1-1,binom_(n2+j,j)*f0(?c,(j+n2,x))
2216 *f1((n1-j),?a)*f2(?b))*force2
2217 +sum_(j,0,n2-1,binom_(n1+j,j)*f0(?c,(j+n1,x))
2218 *f1(?a)*f2((n2-j),?b))*force1
2219 The force operation can be executed directly
2220
2221 The next question is how to program this: recursively or linearly
2222 which would require simulation of a recursion. Recursive is clearest
2223 but we need to pass a number of arguments from the calling routine
2224 to the final routine. This is done with AN.SHvar.
2225
2226 We need space for the accumulation of the combinatoric factors.
2227*/
2228
2229int Shuffle(WORD *from1, WORD *from2, WORD *to)
2230{
2231 GETIDENTITY
2232 WORD *t, *fr, *next1, *next2, na, *fn1, *fn2, *tt;
2233 int i, n, n1, n2, j;
2234 LONG combilast;
2235 SHvariables *SH = &(AN.SHvar);
2236 if ( from1 == SH->stop1 && from2 == SH->stop2 ) {
2237 return(FiniShuffle(to));
2238 }
2239 else if ( from1 == SH->stop1 ) {
2240 i = SH->stop2 - from2; t = to; tt = from2; NCOPY(t,tt,i)
2241 return(FiniShuffle(t));
2242 }
2243 else if ( from2 == SH->stop2 ) {
2244 i = SH->stop1 - from1; t = to; tt = from1; NCOPY(t,tt,i)
2245 return(FiniShuffle(t));
2246 }
2247/*
2248 Compare lead arguments
2249*/
2250 if ( AreArgsEqual(from1,from2) ) {
2251/*
2252 First find out how many of each
2253*/
2254 next1 = from1; n1 = 1; NEXTARG(next1)
2255 while ( ( next1 < SH->stop1 ) && AreArgsEqual(from1,next1) ) {
2256 n1++; NEXTARG(next1)
2257 }
2258 next2 = from2; n2 = 1; NEXTARG(next2)
2259 while ( ( next2 < SH->stop2 ) && AreArgsEqual(from2,next2) ) {
2260 n2++; NEXTARG(next2)
2261 }
2262 combilast = SH->combilast;
2263/*
2264 +binom_(n1+n2,n1)*f0(?c,(n1+n2,x))*f1(?a)*f2(?b)
2265*/
2266 t = to;
2267 n = n1 + n2;
2268 while ( --n >= 0 ) { fr = from1; CopyArg(t,fr) }
2269 if ( GetBinom((UWORD *)(t),&na,n1+n2,n1) ) goto shuffcall;
2270 if ( combilast + AN.SHcombi[combilast] + na + 2 >= AN.SHcombisize ) {
2271/*
2272 We need more memory in this stack. Fortunately this is the
2273 only place where we have to do this, because the other factors
2274 are definitely smaller.
2275 Layout: size, LongInteger, size, LongInteger, .....
2276 We start pointing at the last one.
2277*/
2278 UWORD *combi = (UWORD *)Malloc1(2*AN.SHcombisize*2,"AN.SHcombi");
2279 LONG jj;
2280 for ( jj = 0; jj < AN.SHcombisize; jj++ ) combi[jj] = AN.SHcombi[jj];
2281 AN.SHcombisize *= 2;
2282 M_free(AN.SHcombi,"AN.SHcombi");
2283 AN.SHcombi = combi;
2284 }
2285 if ( MulLong((UWORD *)(AN.SHcombi+combilast+1),AN.SHcombi[combilast],
2286 (UWORD *)(t),na,
2287 (UWORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+2),
2288 (WORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+1)) ) goto shuffcall;
2289 SH->combilast = combilast + AN.SHcombi[combilast] + 1;
2290 if ( next1 >= SH->stop1 ) {
2291 fr = next2; i = SH->stop2 - fr;
2292 NCOPY(t,fr,i)
2293 if ( FiniShuffle(t) ) goto shuffcall;
2294 }
2295 else if ( next2 >= SH->stop2 ) {
2296 fr = next1; i = SH->stop1 - fr;
2297 NCOPY(t,fr,i)
2298 if ( FiniShuffle(t) ) goto shuffcall;
2299 }
2300 else {
2301 if ( Shuffle(next1,next2,t) ) goto shuffcall;
2302 }
2303 SH->combilast = combilast;
2304/*
2305 +sum_(j,0,n1-1,binom_(n2+j,j)*f0(?c,(j+n2,x))
2306 *f1((n1-j),?a)*f2(?b))*force2
2307*/
2308 if ( next2 < SH->stop2 ) {
2309 t = to;
2310 n = n2;
2311 while ( --n >= 0 ) { fr = from1; CopyArg(t,fr) }
2312 for ( j = 0; j < n1; j++ ) {
2313 if ( GetBinom((UWORD *)(t),&na,n2+j,j) ) goto shuffcall;
2314 if ( MulLong((UWORD *)(AN.SHcombi+combilast+1),AN.SHcombi[combilast],
2315 (UWORD *)(t),na,
2316 (UWORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+2),
2317 (WORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+1)) ) goto shuffcall;
2318 SH->combilast = combilast + AN.SHcombi[combilast] + 1;
2319 if ( j > 0 ) { fr = from1; CopyArg(t,fr) }
2320 fn2 = next2; tt = t;
2321 CopyArg(tt,fn2)
2322
2323 if ( fn2 >= SH->stop2 ) {
2324 n = n1-j;
2325 while ( --n >= 0 ) { fr = from1; CopyArg(tt,fr) }
2326 fr = next1; i = SH->stop1 - fr;
2327 NCOPY(tt,fr,i)
2328 if ( FiniShuffle(tt) ) goto shuffcall;
2329 }
2330 else {
2331 n = j; fn1 = from1; while ( --n >= 0 ) { NEXTARG(fn1) }
2332 if ( Shuffle(fn1,fn2,tt) ) goto shuffcall;
2333 }
2334 SH->combilast = combilast;
2335 }
2336 }
2337/*
2338 +sum_(j,0,n2-1,binom_(n1+j,j)*f0(?c,(j+n1,x))
2339 *f1(?a)*f2((n2-j),?b))*force1
2340*/
2341 if ( next1 < SH->stop1 ) {
2342 t = to;
2343 n = n1;
2344 while ( --n >= 0 ) { fr = from1; CopyArg(t,fr) }
2345 for ( j = 0; j < n2; j++ ) {
2346 if ( GetBinom((UWORD *)(t),&na,n1+j,j) ) goto shuffcall;
2347 if ( MulLong((UWORD *)(AN.SHcombi+combilast+1),AN.SHcombi[combilast],
2348 (UWORD *)(t),na,
2349 (UWORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+2),
2350 (WORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+1)) ) goto shuffcall;
2351 SH->combilast = combilast + AN.SHcombi[combilast] + 1;
2352 if ( j > 0 ) { fr = from1; CopyArg(t,fr) }
2353 fn1 = next1; tt = t;
2354 CopyArg(tt,fn1)
2355
2356 if ( fn1 >= SH->stop1 ) {
2357 n = n2-j;
2358 while ( --n >= 0 ) { fr = from1; CopyArg(tt,fr) }
2359 fr = next2; i = SH->stop2 - fr;
2360 NCOPY(tt,fr,i)
2361 if ( FiniShuffle(tt) ) goto shuffcall;
2362 }
2363 else {
2364 n = j; fn2 = from2; while ( --n >= 0 ) { NEXTARG(fn2) }
2365 if ( Shuffle(fn1,fn2,tt) ) goto shuffcall;
2366 }
2367 SH->combilast = combilast;
2368 }
2369 }
2370 }
2371 else {
2372/*
2373 Argument from first list
2374*/
2375 t = to;
2376 fr = from1;
2377 CopyArg(t,fr)
2378 if ( fr >= SH->stop1 ) {
2379 fr = from2; i = SH->stop2 - fr;
2380 NCOPY(t,fr,i)
2381 if ( FiniShuffle(t) ) goto shuffcall;
2382 }
2383 else {
2384 if ( Shuffle(fr,from2,t) ) goto shuffcall;
2385 }
2386/*
2387 Argument from second list
2388*/
2389 t = to;
2390 fr = from2;
2391 CopyArg(t,fr)
2392 if ( fr >= SH->stop2 ) {
2393 fr = from1; i = SH->stop1 - fr;
2394 NCOPY(t,fr,i)
2395 if ( FiniShuffle(t) ) goto shuffcall;
2396 }
2397 else {
2398 if ( Shuffle(from1,fr,t) ) goto shuffcall;
2399 }
2400 }
2401 return(0);
2402shuffcall:
2403 MesCall("Shuffle");
2404 return(-1);
2405}
2406
2407/*
2408 #] Shuffle :
2409 #[ FinishShuffle :
2410
2411 The complications here are:
2412 1: We want to save space. We put the output term in 'out' straight
2413 on top of what we produced thusfar. We have to copy the early
2414 piece because once the term goes back to Generator, Normalize can
2415 change it in situ
2416 2: There can be other occurrence of the function between the two
2417 that we did. For shuffles that isn't likely, but we use this
2418 routine also for the stuffles and there it can happen.
2419*/
2420
2421int FinishShuffle(WORD *fini)
2422{
2423 GETIDENTITY
2424 WORD *t, *t1, *oldworkpointer = AT.WorkPointer, *tcoef, ntcoef, *out;
2425 int i;
2426 SHvariables *SH = &(AN.SHvar);
2427 SH->outfun[1] = fini - SH->outfun;
2428 if ( functions[SH->outfun[0]-FUNCTION].symmetric != 0 )
2429 SH->outfun[2] |= DIRTYSYMFLAG;
2430 out = fini; i = fini - SH->outterm; t = SH->outterm;
2431 NCOPY(fini,t,i)
2432 t = SH->stop1;
2433 t1 = t + t[1];
2434 while ( t1 < SH->stop2 ) { t = t1; t1 = t + t[1]; }
2435 t1 = SH->stop1;
2436 while ( t1 < t ) *fini++ = *t1++;
2437 t = SH->stop2;
2438 while ( t < SH->incoef ) *fini++ = *t++;
2439 tcoef = fini;
2440 ntcoef = SH->nincoef;
2441 i = ABS(ntcoef);
2442 NCOPY(fini,t,i);
2443 ntcoef = REDLENG(ntcoef);
2444 Mully(BHEAD (UWORD *)tcoef,&ntcoef,
2445 (UWORD *)(AN.SHcombi+SH->combilast+1),AN.SHcombi[SH->combilast]);
2446 ntcoef = INCLENG(ntcoef);
2447 fini = tcoef + ABS(ntcoef);
2448 if ( ( ( SH->option & 2 ) != 0 ) && ( ( SH->option & 256 ) != 0 ) ) ntcoef = -ntcoef;
2449 fini[-1] = ntcoef;
2450 i = *out = fini - out;
2451/*
2452 Now check whether we have to do more
2453*/
2454 AT.WorkPointer = out + *out;
2455 if ( ( SH->option & 1 ) == 1 ) {
2456 if ( Generator(BHEAD out,SH->level) ) goto Finicall;
2457 }
2458 else {
2459 if ( DoShtuffle(out,SH->level,SH->thefunction,SH->option) ) goto Finicall;
2460 }
2461 AT.WorkPointer = oldworkpointer;
2462 return(0);
2463Finicall:
2464 AT.WorkPointer = oldworkpointer;
2465 MesCall("FinishShuffle");
2466 return(-1);
2467}
2468
2469/*
2470 #] FinishShuffle :
2471 #[ DoStuffle :
2472
2473 Stuffling is a variation of shuffling.
2474 In the stuffling we insist that the arguments are (short) integers. nonzero.
2475 The stuffle sum is x st y = sig_(x)*sig_(y)*(abs(x)+abs(y))
2476 The way we do this is:
2477 1: count the arguments in each function: n1, n2
2478 2: take the minimum minval = min(n1,n2).
2479 3: for ( j = 0; j <= min; j++ ) take j elements in each of the lists.
2480 4: the j+1 groups of remaining arguments have to each be shuffled
2481 5: the j selected pairs have to be stuffle added.
2482 We can use many of the shuffle things.
2483 Considering the recursive nature of the generation we actually don't
2484 need to know n1, n2, minval.
2485*/
2486
2487WORD DoStuffle(WORD *term, WORD level, WORD fun, WORD option)
2488{
2489 GETIDENTITY
2490 SHvariables SHback, *SH = &(AN.SHvar);
2491 WORD *t1, *t2, *tstop, *t1stop, *t2stop, ncoef, n = fun, *to, *from;
2492 WORD *r1, *r2;
2493 int i, error;
2494 LONG k;
2495 UWORD *newcombi;
2496#ifdef NEWCODE
2497 WORD *rr1, *rr2, i1, i2;
2498#endif
2499 if ( n < 0 ) {
2500 if ( ( n = DolToFunction(BHEAD -n) ) == 0 ) {
2501 MLOCK(ErrorMessageLock);
2502 MesPrint("$-variable in merge statement did not evaluate to a function.");
2503 MUNLOCK(ErrorMessageLock);
2504 return(1);
2505 }
2506 }
2507 if ( AT.WorkPointer + 3*(*term) + AM.MaxTal > AT.WorkTop ) {
2508 MLOCK(ErrorMessageLock);
2509 MesWork();
2510 MUNLOCK(ErrorMessageLock);
2511 return(-1);
2512 }
2513
2514 tstop = term + *term;
2515 ncoef = tstop[-1];
2516 tstop -= ABS(ncoef);
2517 t1 = term + 1;
2518retry1:;
2519 while ( t1 < tstop ) {
2520 if ( ( *t1 == n ) && ( t1+t1[1] < tstop ) && ( t1[1] > FUNHEAD ) ) {
2521 t2 = t1 + t1[1];
2522 if ( t2 >= tstop ) {
2523 return(Generator(BHEAD term,level));
2524 }
2525retry2:;
2526 while ( t2 < tstop ) {
2527 if ( ( *t2 == n ) && ( t2[1] > FUNHEAD ) ) break;
2528 t2 += t2[1];
2529 }
2530 if ( t2 < tstop ) break;
2531 }
2532 t1 += t1[1];
2533 }
2534 if ( t1 >= tstop ) {
2535 return(Generator(BHEAD term,level));
2536 }
2537/*
2538 Next we have to check that the arguments are of the correct type
2539 At the same time we can count them.
2540*/
2541#ifndef NEWCODE
2542 t1stop = t1 + t1[1];
2543 r1 = t1 + FUNHEAD;
2544 while ( r1 < t1stop ) {
2545 if ( *r1 != -SNUMBER ) break;
2546 if ( r1[1] == 0 ) break;
2547 r1 += 2;
2548 }
2549 if ( r1 < t1stop ) { t1 = t2; goto retry1; }
2550 t2stop = t2 + t2[1];
2551 r2 = t2 + FUNHEAD;
2552 while ( r2 < t2stop ) {
2553 if ( *r2 != -SNUMBER ) break;
2554 if ( r2[1] == 0 ) break;
2555 r2 += 2;
2556 }
2557 if ( r2 < t2stop ) { t2 = t2 + t2[1]; goto retry2; }
2558#else
2559 t1stop = t1 + t1[1];
2560 r1 = t1 + FUNHEAD;
2561 while ( r1 < t1stop ) {
2562 if ( *r1 == -SNUMBER ) {
2563 if ( r1[1] == 0 ) break;
2564 r1 += 2; continue;
2565 }
2566 else if ( *r1 == -SYMBOL ) {
2567 if ( ( symbols[r1[1]].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
2568 break;
2569 r1 += 2; continue;
2570 }
2571 if ( *r1 > 0 && *r1 == r1[ARGHEAD]+ARGHEAD ) {
2572 if ( ABS(r1[r1[0]-1]) == r1[0]-ARGHEAD-1 ) {}
2573 else if ( r1[ARGHEAD+1] == SYMBOL ) {
2574 rr1 = r1 + ARGHEAD + 3;
2575 i1 = rr1[-1]-2;
2576 while ( i1 > 0 ) {
2577 if ( ( symbols[*rr1].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
2578 break;
2579 i1 -= 2; rr1 += 2;
2580 }
2581 if ( i1 > 0 ) break;
2582 }
2583 else break;
2584 rr1 = r1+*r1-1;
2585 i1 = (ABS(*rr1)-1)/2;
2586 while ( i1 > 1 ) {
2587 if ( rr1[-1] ) break;
2588 i1--; rr1--;
2589 }
2590 if ( i1 > 1 || rr1[-1] != 1 ) break;
2591 r1 += *r1;
2592 }
2593 else break;
2594 }
2595 if ( r1 < t1stop ) { t1 = t2; goto retry1; }
2596 t2stop = t2 + t2[1];
2597 r2 = t2 + FUNHEAD;
2598
2599 while ( r2 < t2stop ) {
2600 if ( *r2 == -SNUMBER ) {
2601 if ( r2[1] == 0 ) break;
2602 r2 += 2; continue;
2603 }
2604 else if ( *r2 == -SYMBOL ) {
2605 if ( ( symbols[r2[1]].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
2606 break;
2607 r2 += 2; continue;
2608 }
2609 if ( *r2 > 0 && *r2 == r2[ARGHEAD]+ARGHEAD ) {
2610 if ( ABS(r2[r2[0]-1]) == r2[0]-ARGHEAD-1 ) {}
2611 else if ( r2[ARGHEAD+1] == SYMBOL ) {
2612 rr2 = r2 + ARGHEAD + 3;
2613 i2 = rr2[-1]-2;
2614 while ( i2 > 0 ) {
2615 if ( ( symbols[*rr2].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
2616 break;
2617 i2 -= 2; rr2 += 2;
2618 }
2619 if ( i2 > 0 ) break;
2620 }
2621 else break;
2622 rr2 = r2+*r2-1;
2623 i2 = (ABS(*rr2)-1)/2;
2624 while ( i2 > 1 ) {
2625 if ( rr2[-1] ) break;
2626 i2--; rr2--;
2627 }
2628 if ( i2 > 1 || rr2[-1] != 1 ) break;
2629 r2 += *r2;
2630 }
2631 else break;
2632 }
2633 if ( r2 < t2stop ) { t2 = t2 + t2[1]; goto retry2; }
2634#endif
2635/*
2636 OK, now we got two objects that can be used.
2637*/
2638 *AN.RepPoint = 1;
2639
2640 SHback = AN.SHvar;
2641 SH->finishuf = &FinishStuffle;
2642 SH->do_uffle = &DoStuffle;
2643 SH->outterm = AT.WorkPointer;
2644 AT.WorkPointer += *term;
2645 SH->ststop1 = t1 + t1[1];
2646 SH->ststop2 = t2 + t2[1];
2647 SH->thefunction = n;
2648 SH->option = option;
2649 SH->level = level;
2650 SH->incoef = tstop;
2651 SH->nincoef = ncoef;
2652 if ( AN.SHcombi == 0 || AN.SHcombisize == 0 ) {
2653 AN.SHcombisize = 200;
2654 AN.SHcombi = (UWORD *)Malloc1(AN.SHcombisize*sizeof(UWORD),"AN.SHcombi");
2655 SH->combilast = 0;
2656 SHback.combilast = 0;
2657 }
2658 else {
2659 SH->combilast += AN.SHcombi[SH->combilast]+1;
2660 if ( SH->combilast >= AN.SHcombisize - 100 ) {
2661 newcombi = (UWORD *)Malloc1(2*AN.SHcombisize*sizeof(UWORD),"AN.SHcombi");
2662 for ( k = 0; k < AN.SHcombisize; k++ ) newcombi[k] = AN.SHcombi[k];
2663 M_free(AN.SHcombi,"AN.SHcombi");
2664 AN.SHcombi = newcombi;
2665 AN.SHcombisize *= 2;
2666 }
2667 }
2668 AN.SHcombi[SH->combilast] = 1;
2669 AN.SHcombi[SH->combilast+1] = 1;
2670
2671 i = t1-term; to = SH->outterm; from = term;
2672 NCOPY(to,from,i)
2673 SH->outfun = to;
2674 for ( i = 0; i < FUNHEAD; i++ ) { *to++ = t1[i]; }
2675
2676 error = Stuffle(t1+FUNHEAD,t2+FUNHEAD,to);
2677
2678 AT.WorkPointer = SH->outterm;
2679 AN.SHvar = SHback;
2680 if ( error ) {
2681 MesCall("DoStuffle");
2682 return(-1);
2683 }
2684 return(0);
2685}
2686
2687/*
2688 #] DoStuffle :
2689 #[ Stuffle :
2690
2691 The way to generate the stuffles
2692 1: select an argument in the first list (for(j1=0;j1<last;j1++))
2693 2: select an argument in the second list (for(j2=0;j2<last;j2++))
2694 3: put values for SH->ststop1 and SH->ststop2 at these arguments.
2695 4: generate all shuffles of the arguments in front.
2696 5: Then put the stuffle sum of arg(j1) and arg(j2)
2697 6: Then continue calling Stuffle
2698 7: Once one gets exhausted, we can clean up the list and call FinishShuffle
2699 8: if ( ( SH->option & 2 ) != 0 ) the stuffle sum is negative.
2700*/
2701
2702int Stuffle(WORD *from1, WORD *from2, WORD *to)
2703{
2704 GETIDENTITY
2705 WORD *t, *tf, *next1, *next2, *st1, *st2, *save1, *save2;
2706 SHvariables *SH = &(AN.SHvar);
2707 int i, retval;
2708/*
2709 First the special cases (exhausted list(s)):
2710*/
2711 save1 = SH->stop1; save2 = SH->stop2;
2712 if ( from1 >= SH->ststop1 && from2 == SH->ststop2 ) {
2713 SH->stop1 = SH->ststop1;
2714 SH->stop2 = SH->ststop2;
2715 retval = FinishShuffle(to);
2716 SH->stop1 = save1; SH->stop2 = save2;
2717 return(retval);
2718 }
2719 else if ( from1 >= SH->ststop1 ) {
2720 i = SH->ststop2 - from2; t = to; tf = from2; NCOPY(t,tf,i)
2721 SH->stop1 = SH->ststop1;
2722 SH->stop2 = SH->ststop2;
2723 retval = FinishShuffle(t);
2724 SH->stop1 = save1; SH->stop2 = save2;
2725 return(retval);
2726 }
2727 else if ( from2 >= SH->ststop2 ) {
2728 i = SH->ststop1 - from1; t = to; tf = from1; NCOPY(t,tf,i)
2729 SH->stop1 = SH->ststop1;
2730 SH->stop2 = SH->ststop2;
2731 retval = FinishShuffle(t);
2732 SH->stop1 = save1; SH->stop2 = save2;
2733 return(retval);
2734 }
2735/*
2736 Now the case that we have no stuffle sums.
2737*/
2738 SH->stop1 = SH->ststop1;
2739 SH->stop2 = SH->ststop2;
2740 SH->finishuf = &FinishShuffle;
2741 if ( Shuffle(from1,from2,to) ) goto stuffcall;
2742 SH->finishuf = &FinishStuffle;
2743/*
2744 Now we have to select a pair, one from 1 and one from 2.
2745*/
2746#ifndef NEWCODE
2747 st1 = from1; next1 = st1+2; /* <----- */
2748#else
2749 st1 = next1 = from1;
2750 NEXTARG(next1)
2751#endif
2752 while ( next1 <= SH->ststop1 ) {
2753#ifndef NEWCODE
2754 st2 = from2; next2 = st2+2; /* <----- */
2755#else
2756 next2 = st2 = from2;
2757 NEXTARG(next2)
2758#endif
2759 while ( next2 <= SH->ststop2 ) {
2760 SH->stop1 = st1;
2761 SH->stop2 = st2;
2762 if ( st1 == from1 && st2 == from2 ) {
2763 t = to;
2764#ifndef NEWCODE
2765 *t++ = -SNUMBER; *t++ = StuffAdd(st1[1],st2[1]);
2766#else
2767 t = StuffRootAdd(st1,st2,t);
2768#endif
2769 SH->option ^= 256;
2770 if ( Stuffle(next1,next2,t) ) goto stuffcall;
2771 SH->option ^= 256;
2772 }
2773 else if ( st1 == from1 ) {
2774 i = st2-from2;
2775 t = to; tf = from2; NCOPY(t,tf,i)
2776#ifndef NEWCODE
2777 *t++ = -SNUMBER; *t++ = StuffAdd(st1[1],st2[1]);
2778#else
2779 t = StuffRootAdd(st1,st2,t);
2780#endif
2781 SH->option ^= 256;
2782 if ( Stuffle(next1,next2,t) ) goto stuffcall;
2783 SH->option ^= 256;
2784 }
2785 else if ( st2 == from2 ) {
2786 i = st1-from1;
2787 t = to; tf = from1; NCOPY(t,tf,i)
2788#ifndef NEWCODE
2789 *t++ = -SNUMBER; *t++ = StuffAdd(st1[1],st2[1]);
2790#else
2791 t = StuffRootAdd(st1,st2,t);
2792#endif
2793 SH->option ^= 256;
2794 if ( Stuffle(next1,next2,t) ) goto stuffcall;
2795 SH->option ^= 256;
2796 }
2797 else {
2798 if ( Shuffle(from1,from2,to) ) goto stuffcall;
2799 }
2800#ifndef NEWCODE
2801 st2 = next2; next2 += 2; /* <----- */
2802#else
2803 st2 = next2;
2804 NEXTARG(next2)
2805#endif
2806 }
2807#ifndef NEWCODE
2808 st1 = next1; next1 += 2; /* <----- */
2809#else
2810 st1 = next1;
2811 NEXTARG(next1)
2812#endif
2813 }
2814 SH->stop1 = save1; SH->stop2 = save2;
2815 return(0);
2816stuffcall:;
2817 MesCall("Stuffle");
2818 return(-1);
2819}
2820
2821/*
2822 #] Stuffle :
2823 #[ FinishStuffle :
2824
2825 The program only comes here from the Shuffle routine.
2826 It should add the stuffle sum and then call Stuffle again.
2827*/
2828
2829int FinishStuffle(WORD *fini)
2830{
2831 GETIDENTITY
2832 SHvariables *SH = &(AN.SHvar);
2833#ifdef NEWCODE
2834 WORD *next1 = SH->stop1, *next2 = SH->stop2;
2835 fini = StuffRootAdd(next1,next2,fini);
2836#else
2837 *fini++ = -SNUMBER; *fini++ = StuffAdd(SH->stop1[1],SH->stop2[1]);
2838#endif
2839 SH->option ^= 256;
2840#ifdef NEWCODE
2841 NEXTARG(next1)
2842 NEXTARG(next2)
2843 if ( Stuffle(next1,next2,fini) ) goto stuffcall;
2844#else
2845 if ( Stuffle(SH->stop1+2,SH->stop2+2,fini) ) goto stuffcall;
2846#endif
2847 SH->option ^= 256;
2848 return(0);
2849stuffcall:;
2850 MesCall("FinishStuffle");
2851 return(-1);
2852}
2853
2854/*
2855 #] FinishStuffle :
2856 #[ StuffRootAdd :
2857
2858 Makes the stuffle sum of two arguments.
2859 The arguments can be of one of three types:
2860 1: -SNUMBER,num
2861 2: -SYMBOL,symbol
2862 3: Numerical (long) argument.
2863 4: Generic argument with (only) symbols that are roots of unity and
2864 a coefficient.
2865 We have excluded the case that both t1 and t2 are of type 1:
2866 The output should be written to 'to' and the new fill position should
2867 be the return value.
2868 `to' is inside the workspace.
2869
2870 The stuffle sum is sig_(t2)*t1+sig_(t1)*t2
2871 or sig_(t1)*sig_(t2)*(abs_(t1)+abs_(t2))
2872*/
2873
2874#ifdef NEWCODE
2875
2876WORD *StuffRootAdd(WORD *t1, WORD *t2, WORD *to)
2877{
2878 int type1, type2, type3, sgn, sgn1, sgn2, sgn3, pow, root, nosymbols, i;
2879 WORD *tt1, *tt2, it1, it2, *t3, *r, size1, size2, size3;
2880 WORD scratch[2];
2881 LONG x;
2882 if ( *t1 == -SNUMBER ) { type1 = 1; if ( t1[1] < 0 ) sgn1 = -1; else sgn1 = 1; }
2883 else if ( *t1 == -SYMBOL ) { type1 = 2; sgn1 = 1; }
2884 else if ( ABS(t1[*t1-1]) == *t1-ARGHEAD-1 ) {
2885 type1 = 3; if ( t1[*t1-1] < 0 ) sgn1 = -1; else sgn1 = 1; }
2886 else { type1 = 4; if ( t1[*t1-1] < 0 ) sgn1 = -1; else sgn1 = 1; }
2887 if ( *t2 == -SNUMBER ) { type2 = 1; if ( t2[1] < 0 ) sgn2 = -1; else sgn2 = 1; }
2888 else if ( *t2 == -SYMBOL ) { type2 = 2; sgn2 = 1; }
2889 else if ( ABS(t2[*t2-1]) == *t2-ARGHEAD-1 ) {
2890 type2 = 3; if ( t2[*t2-1] < 0 ) sgn2 = -1; else sgn2 = 1; }
2891 else { type2 = 4; if ( t2[*t2-1] < 0 ) sgn2 = -1; else sgn2 = 1; }
2892 if ( type1 > type2 ) {
2893 t3 = t1; t1 = t2; t2 = t3;
2894 type3 = type1; type1 = type2; type2 = type3;
2895 sgn3 = sgn1; sgn1 = sgn2; sgn2 = sgn3;
2896 }
2897 nosymbols = 1; sgn3 = 1;
2898 switch ( type1 ) {
2899 case 1:
2900 if ( type2 == 1 ) {
2901 x = sgn2 * t1[1];
2902 x += sgn1 * t2[1];
2903 if ( x > MAXPOSITIVE || x < -(MAXPOSITIVE+1) ) {
2904 if ( x < 0 ) { sgn1 = -3; x = -x; }
2905 else sgn1 = 3;
2906 *to++ = ARGHEAD+4;
2907 *to++ = 0;
2908 FILLARG(to)
2909 *to++ = 4; *to++ = (UWORD)x; *to++ = 1; *to++ = sgn1;
2910 }
2911 else { *to++ = -SNUMBER; *to++ = (WORD)x; }
2912 }
2913 else if ( type2 == 2 ) {
2914 *to++ = ARGHEAD+8; *to++ = 0; FILLARG(to)
2915 *to++ = 8; *to++ = SYMBOL; *to++ = 4; *to++ = t2[1]; *to++ = 1;
2916 *to++ = ABS(t1[1])+1;
2917 *to++ = 1;
2918 *to++ = 3*sgn1;
2919 }
2920 else if ( type2 == 3 ) {
2921 tt1 = (WORD *)scratch; tt1[0] = ABS(t1[1]); size1 = 1;
2922 tt2 = t2+ARGHEAD+1; size2 = (ABS(t2[*t2-1])-1)/2;
2923 t3 = to;
2924 *to++ = 0; *to++ = 0; FILLARG(to) *to++ = 0;
2925 goto DoCoeffi;
2926 }
2927 else {
2928/*
2929 t1 is (short) numeric, t2 has the symbol(s).
2930*/
2931 tt1 = (WORD *)scratch; tt1[0] = ABS(t1[1]); size1 = 1;
2932 tt2 = t2+ARGHEAD+1; tt2 += tt2[1]; size2 = (ABS(t2[*t2-1])-1)/2;
2933 t3 = to; i = tt2 - t2; r = t2;
2934 NCOPY(to,r,i)
2935 nosymbols = 0;
2936 goto DoCoeffi;
2937 }
2938 break;
2939 case 2:
2940 if ( type2 == 2 ) {
2941 if ( t1[1] == t2[1] ) {
2942 if ( ( symbols[t1[1]].maxpower == 4 )
2943 && ( ( symbols[t1[1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) ) {
2944 *to++ = -SNUMBER; *to++ = -2;
2945 }
2946 else if ( symbols[t1[1]].maxpower == 2 ) {
2947 *to++ = -SNUMBER; *to++ = 2;
2948 }
2949 else {
2950 *to++ = ARGHEAD+8; *to++ = 0; FILLARG(to)
2951 *to++ = 8; *to++ = SYMBOL; *to++ = 4;
2952 *to++ = t1[1]; *to++ = 2;
2953 *to++ = 2; *to++ = 1; *to++ = 3;
2954 }
2955 }
2956 else {
2957 *to++ = ARGHEAD+10; *to++ = 0; FILLARG(to)
2958 *to++ = 10; *to++ = SYMBOL; *to++ = 6;
2959 if ( t1[1] < t2[1] ) {
2960 *to++ = t1[1]; *to++ = 1; *to++ = t2[1]; *to++ = 1;
2961 }
2962 else {
2963 *to++ = t2[1]; *to++ = 1; *to++ = t1[1]; *to++ = 1;
2964 }
2965 *to++ = 2; *to++ = 1; *to++ = 3;
2966 }
2967 }
2968 else if ( type2 == 3 ) {
2969 t3 = to;
2970 *to++ = 0; *to++ = 0; FILLARG(to) *to++ = 0;
2971 *to++ = SYMBOL; *to++ = 4; *to++ = t1[1]; *to++ = 1;
2972 tt1 = scratch; tt1[1] = 1; size1 = 1;
2973 tt2 = t2+ARGHEAD+1; size2 = (ABS(t2[*t2-1])-1)/2;
2974 nosymbols = 0;
2975 goto DoCoeffi;
2976 }
2977 else {
2978 tt1 = scratch; tt1[0] = 1; size1 = 1;
2979 t3 = to;
2980 *to++ = 0; *to++ = 0; FILLARG(to) *to++ = 0;
2981 *to++ = SYMBOL; *to++ = 0;
2982 tt2 = t2 + ARGHEAD+3; it2 = tt2[-1]-2;
2983 while ( it2 > 0 ) {
2984 if ( *tt2 == t1[1] ) {
2985 pow = tt2[1]+1;
2986 root = symbols[*tt2].maxpower;
2987 if ( pow >= root ) pow -= root;
2988 if ( ( symbols[*tt2].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
2989 if ( ( root & 1 ) == 0 && pow >= root/2 ) {
2990 pow -= root/2; sgn3 = -sgn3;
2991 }
2992 }
2993 if ( pow != 0 ) {
2994 *to++ = *tt2; *to++ = pow;
2995 }
2996 tt2 += 2; it2 -= 2;
2997 break;
2998 }
2999 else if ( t1[1] < *tt2 ) {
3000 *to++ = t1[1]; *to++ = 1; break;
3001 }
3002 else {
3003 *to++ = *tt2++; *to++ = *tt2++; it2 -= 2;
3004 if ( it2 <= 0 ) { *to++ = t1[1]; *to++ = 1; }
3005 }
3006 }
3007 while ( it2 > 0 ) { *to++ = *tt2++; *to++ = *tt2++; it2 -= 2; }
3008 if ( (to - t3) > ARGHEAD+3 ) {
3009 t3[ARGHEAD+2] = (to-t3)-ARGHEAD-1; /* size of the SYMBOL field */
3010 nosymbols = 0;
3011 }
3012 else {
3013 to = t3+ARGHEAD+1; /* no SYMBOL field */
3014 }
3015 size2 = (ABS(t2[*t2-1])-1)/2;
3016 goto DoCoeffi;
3017 }
3018 break;
3019 case 3:
3020 if ( type2 == 3 ) {
3021/*
3022 Both are numeric
3023*/
3024 tt1 = t1+ARGHEAD+1; size1 = (ABS(t1[*t1-1])-1)/2;
3025 tt2 = t2+ARGHEAD+1; size2 = (ABS(t2[*t2-1])-1)/2;
3026 t3 = to;
3027 *to++ = 0; *to++ = 0; FILLARG(to) *to++ = 0;
3028 goto DoCoeffi;
3029 }
3030 else {
3031/*
3032 t1 is (long) numeric, t2 has the symbol(s).
3033*/
3034 tt1 = t1+ARGHEAD+1; size1 = (ABS(t1[*t1-1])-1)/2;
3035 tt2 = t2+ARGHEAD+1; tt2 += tt2[1]; size2 = (ABS(t2[*t2-1])-1)/2;
3036 t3 = to; i = tt2 - t2; r = t2;
3037 NCOPY(to,r,i)
3038 nosymbols = 0;
3039 goto DoCoeffi;
3040 }
3041 break;
3042 case 4:
3043/*
3044 Both have roots of unity
3045 1: Merge the lists and simplify if possible
3046*/
3047 tt1 = t1+ARGHEAD+3; it1 = tt1[-1]-2;
3048 tt2 = t2+ARGHEAD+3; it2 = tt2[-1]-2;
3049 t3 = to;
3050 *to++ = 0; *to++ = 0; FILLARG(to)
3051 *to++ = 0; *to++ = SYMBOL; *to++ = 0;
3052 while ( it1 > 0 && it2 > 0 ) {
3053 if ( *tt1 == *tt2 ) {
3054 pow = tt1[1]+tt2[1];
3055 root = symbols[*tt1].maxpower;
3056 if ( pow >= root ) pow -= root;
3057 if ( ( symbols[*tt1].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
3058 if ( ( root & 1 ) == 0 && pow >= root/2 ) {
3059 pow -= root/2; sgn3 = -sgn3;
3060 }
3061 }
3062 if ( pow != 0 ) {
3063 *to++ = *tt1; *to++ = pow;
3064 }
3065 tt1 += 2; tt2 += 2; it1 -= 2; it2 -= 2;
3066 }
3067 else if ( *tt1 < *tt2 ) {
3068 *to++ = *tt1++; *to++ = *tt1++; it1 -= 2;
3069 }
3070 else {
3071 *to++ = *tt2++; *to++ = *tt2++; it2 -= 2;
3072 }
3073 }
3074 while ( it1 > 0 ) { *to++ = *tt1++; *to++ = *tt1++; it1 -= 2; }
3075 while ( it2 > 0 ) { *to++ = *tt2++; *to++ = *tt2++; it2 -= 2; }
3076 if ( (to - t3) > ARGHEAD+3 ) {
3077 t3[ARGHEAD+2] = (to-t3)-ARGHEAD-1; /* size of the SYMBOL field */
3078 nosymbols = 0;
3079 }
3080 else {
3081 to = t3+ARGHEAD+1; /* no SYMBOL field */
3082 }
3083 size1 = (ABS(t1[*t1-1])-1)/2;
3084 size2 = (ABS(t2[*t2-1])-1)/2;
3085/*
3086 Now tt1 and tt2 are pointing at their coefficients.
3087 sgn1 is the sign of 1, sgn2 is the sign of 2 and sgn3 is an extra
3088 overall sign.
3089*/
3090DoCoeffi:
3091 if ( AddLong((UWORD *)tt1,size1,(UWORD *)tt2,size2,(UWORD *)to,&size3) ) {
3092 MLOCK(ErrorMessageLock);
3093 MesPrint("Called from StuffRootAdd");
3094 MUNLOCK(ErrorMessageLock);
3095 Terminate(-1);
3096 }
3097 sgn = sgn1*sgn2*sgn3;
3098 if ( nosymbols && size3 == 1 ) {
3099 if ( (UWORD)(to[0]) <= MAXPOSITIVE && sgn > 0 ) {
3100 sgn1 = to[0];
3101 to = t3; *to++ = -SNUMBER; *to++ = sgn1;
3102 }
3103 else if ( (UWORD)(to[0]) <= (MAXPOSITIVE+1) && sgn < 0 ) {
3104 sgn1 = to[0];
3105 to = t3; *to++ = -SNUMBER; *to++ = -sgn1;
3106 }
3107 else goto genericcoef;
3108 }
3109 else {
3110genericcoef:
3111 to += size3;
3112 sgn = sgn*(2*size3+1);
3113 *to++ = 1;
3114 while ( size3 > 1 ) { *to++ = 0; size3--; }
3115 *to++ = sgn;
3116 t3[0] = to - t3;
3117 t3[ARGHEAD] = t3[0] - ARGHEAD;
3118 }
3119 break;
3120 }
3121 return(to);
3122}
3123
3124#endif
3125
3126/*
3127 #] StuffRootAdd :
3128*/
WORD CompCoef(WORD *, WORD *)
Definition reken.c:3037
WORD Generator(PHEAD WORD *, WORD)
Definition proces.c:3101
struct PaRtI PARTI
struct PeRmUtEp PERMP