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