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