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