FORM  4.1
reshuf.c
Go to the documentation of this file.
1 
9 /* #[ License : */
10 /*
11  * Copyright (C) 1984-2013 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 ) return(Generator(BHEAD r,level));
1092  if ( i == 0 && CompCoef(term,r) != 0 ) { return(0); }
1093  }
1094  AT.WorkPointer = r;
1095  return(Generator(BHEAD term,level));
1096 }
1097 
1098 /*
1099  #] TryDo :
1100  #] Try Term(s) :
1101  #[ Distribute :
1102  #[ DoDistrib :
1103 
1104  The routine that generates the terms ordered by a distrib_ function.
1105  The presence of a replaceable distrib_ function has been sensed
1106  in the routine TestSub and has been passed on to Generator.
1107  It is then Generator that calls this function in a way that is
1108  similar to calling the trace routines, except for that for the
1109  trace routines and the Levi-Civita tensors the arguments are put
1110  in temporary storage and here we leave them inside the term,
1111  because there is no knowing how long the field will be.
1112 */
1113 
1114 WORD DoDistrib(PHEAD WORD *term, WORD level)
1115 {
1116  GETBIDENTITY
1117  WORD *t, *m, *r = 0, *stop, *tstop, *termout, *endhead, *starttail, *parms;
1118  WORD i, j, k, n, nn, ntype, fun1 = 0, fun2 = 0, typ1 = 0, typ2 = 0;
1119  WORD *arg, *oldwork, *mf, ktype = 0, atype = 0;
1120  WORD sgn, dirtyflag;
1121  AN.TeInFun = AR.TePos = 0;
1122  t = term;
1123  tstop = t + *t;
1124  stop = tstop - ABS(tstop[-1]);
1125  t++;
1126  while ( t < stop ) {
1127  r = t + t[1];
1128  if ( *t == DISTRIBUTION && t[FUNHEAD] == -SNUMBER
1129  && t[FUNHEAD+1] >= -2 && t[FUNHEAD+1] <= 2
1130  && t[FUNHEAD+2] == -SNUMBER
1131  && t[FUNHEAD+4] <= -FUNCTION
1132  && t[FUNHEAD+5] <= -FUNCTION ) {
1133  fun1 = -t[FUNHEAD+4];
1134  fun2 = -t[FUNHEAD+5];
1135  typ1 = functions[fun1-FUNCTION].spec;
1136  typ2 = functions[fun2-FUNCTION].spec;
1137  if ( typ1 > 0 || typ2 > 0 ) {
1138  m = t + FUNHEAD+6;
1139  r = t + t[1];
1140  while ( m < r ) {
1141  if ( *m != -INDEX && *m != -VECTOR && *m != -MINVECTOR )
1142  break;
1143  m += 2;
1144  }
1145  if ( m < r ) {
1146  MLOCK(ErrorMessageLock);
1147  MesPrint("Incompatible function types and arguments in distrib_");
1148  MUNLOCK(ErrorMessageLock);
1149  SETERROR(-1)
1150  }
1151  }
1152  break;
1153  }
1154  t = r;
1155  }
1156  dirtyflag = t[2];
1157  ntype = t[FUNHEAD+1];
1158  n = t[FUNHEAD+3];
1159 /*
1160  t points at the distrib_ function to be expanded.
1161  fun1,fun2 and typ1,typ2 are the two functions and their types.
1162  ntype indicates the action:
1163  0: Make all possible divisions: 2^nargs
1164  1: fun1 should get n arguments: nargs! / ( n! (nargs-n)! )
1165  2: fun2 should get n arguments: nargs! / ( n! (nargs-n)! )
1166  The distiction between 1 and two is for noncommuting objects.
1167  3: fun1 should get n arguments. Super symmetric option.
1168  4: fun2 idem
1169  The super symmetric option involves:
1170  a: arguments get sorted
1171  b: identical arguments are seen as such. Hence not all their
1172  distributions are taken into account. It is as if after the
1173  distrib there is a symmetrize fun1; symmetrize fun2;
1174  c: Hence if the occurrence of each argument is a,b,c,...
1175  and their occurrence in fun1 is a1,b1,c1,... and in fun2
1176  is a2,b2,c2,... then each term is generated (a1+a2)!/a1!/a2!
1177  (b1+b2)!/b1!/b2! (c1+c2)!/c1!/c2! ... times.
1178  d: We have to make an array of occurrences and positions.
1179  e: Then we sort the arguments indirectly.
1180  f: Next we generate the argument lists in the same way as we
1181  generate powers of expressions with binomials. Hence we need
1182  a third array to keep track of the `powers'
1183 */
1184  endhead = t;
1185  starttail = r;
1186  parms = m = t + FUNHEAD+6;
1187  i = 0;
1188  while ( m < r ) { /* Count arguments */
1189  i++;
1190  NEXTARG(m);
1191  }
1192  oldwork = AT.WorkPointer;
1193  arg = AT.WorkPointer + 1;
1194  arg[-1] = 0;
1195  termout = arg + i;
1196  switch ( ntype ) {
1197  case 0: ktype = 1; atype = n < 0 ? 1: 0; n = 0; break;
1198  case 1: ktype = 1; atype = 0; break;
1199  case 2: ktype = 0; atype = 0; break;
1200  case -1: ktype = 1; atype = 1; break;
1201  case -2: ktype = 0; atype = 1; break;
1202  }
1203  do {
1204 /*
1205  All distributions with n elements. We generate the array arg with
1206  all possible 1 and 0 patterns. 1 means in fun1 and 0 means in fun2.
1207 */
1208  if ( n > i ) return(0); /* 0 elements */
1209 
1210  for ( j = 0; j < n; j++ ) arg[j] = 1;
1211  for ( j = n; j < i; j++ ) arg[j] = 0;
1212  for(;;) {
1213  sgn = 0;
1214  t = term;
1215  m = termout;
1216  while ( t < endhead ) *m++ = *t++;
1217  mf = m;
1218  *m++ = fun1;
1219  *m++ = FUNHEAD;
1220  *m++ = dirtyflag;
1221 #if FUNHEAD > 3
1222  k = FUNHEAD -3;
1223  while ( k-- > 0 ) *m++ = 0;
1224 #endif
1225  r = parms;
1226  for ( k = 0; k < i; k++ ) {
1227  if ( arg[k] == ktype ) {
1228  if ( *r <= -FUNCTION ) *m++ = *r++;
1229  else if ( *r < 0 ) {
1230  if ( typ1 > 0 ) {
1231  if ( *r == -MINVECTOR ) sgn ^= 1;
1232  r++;
1233  *m++ = *r++;
1234  }
1235  else { *m++ = *r++; *m++ = *r++; }
1236  }
1237  else {
1238  nn = *r;
1239  NCOPY(m,r,nn);
1240  }
1241  }
1242  else { NEXTARG(r) }
1243  }
1244  mf[1] = WORDDIF(m,mf);
1245  mf = m;
1246  *m++ = fun2;
1247  *m++ = FUNHEAD;
1248  *m++ = dirtyflag;
1249 #if FUNHEAD > 3
1250  k = FUNHEAD -3;
1251  while ( k-- > 0 ) *m++ = 0;
1252 #endif
1253  r = parms;
1254  for ( k = 0; k < i; k++ ) {
1255  if ( arg[k] != ktype ) {
1256  if ( *r <= -FUNCTION ) *m++ = *r++;
1257  else if ( *r < 0 ) {
1258  if ( typ2 > 0 ) {
1259  if ( *r == -MINVECTOR ) sgn ^= 1;
1260  r++;
1261  *m++ = *r++;
1262  }
1263  else { *m++ = *r++; *m++ = *r++; }
1264  }
1265  else {
1266  nn = *r;
1267  NCOPY(m,r,nn);
1268  }
1269  }
1270  else { NEXTARG(r) }
1271  }
1272  mf[1] = WORDDIF(m,mf);
1273 #ifndef NUOVO
1274  if ( atype == 0 ) {
1275  WORD k1,k2;
1276  for ( k = 0; k < i-1; k++ ) {
1277  if ( arg[k] == 0 ) continue;
1278  k1 = 1; k2 = k;
1279  while ( k < i-1 && EqualArg(parms,k,k+1) ) { k++; k1++; }
1280  while ( k2 <= k && arg[k2] == 1 ) k2++;
1281  k2 = k-k2+1;
1282 /*
1283  Now we need k1!/(k2! (k1-k2)!)
1284 */
1285  if ( k2 != k1 && k2 != 0 ) {
1286  if ( GetBinom((UWORD *)m+3,m+2,k1,k2) ) {
1287  MLOCK(ErrorMessageLock);
1288  MesCall("DoDistrib");
1289  MUNLOCK(ErrorMessageLock);
1290  SETERROR(-1)
1291  }
1292  m[1] = ( m[2] < 0 ? -m[2]: m[2] ) + 3;
1293  *m = LNUMBER;
1294  m += m[1];
1295  }
1296  }
1297  }
1298 #endif
1299  r = starttail;
1300  while ( r < tstop ) *m++ = *r++;
1301 
1302  if ( atype ) { /* antisymmetric field */
1303  k = n;
1304  nn = 0;
1305  for ( j = 0; j < i && k > 0; j++ ) {
1306  if ( arg[j] == 1 ) k--;
1307  else nn += k;
1308  }
1309  sgn ^= nn & 1;
1310  }
1311 
1312  if ( sgn ) m[-1] = -m[-1];
1313  *termout = WORDDIF(m,termout);
1314  AT.WorkPointer = m;
1315  if ( AT.WorkPointer > AT.WorkTop ) {
1316  MLOCK(ErrorMessageLock);
1317  MesWork();
1318  MUNLOCK(ErrorMessageLock);
1319  return(-1);
1320  }
1321  *AN.RepPoint = 1;
1322  AR.expchanged = 1;
1323  if ( Generator(BHEAD termout,level) ) Terminate(-1);
1324 #ifndef NUOVO
1325  {
1326  WORD k1;
1327  j = i - 1;
1328  k = 0;
1329 redok: while ( arg[j] == 1 && j >= 0 ) { j--; k++; }
1330  while ( arg[j] == 0 && j >= 0 ) j--;
1331  if ( j < 0 ) break;
1332  k1 = j;
1333  arg[j] = 0;
1334  while ( !atype && EqualArg(parms,j,j+1) ) {
1335  j++;
1336  if ( j >= i - k - 1 ) { j = k1; k++; goto redok; }
1337  arg[j] = 0;
1338  }
1339  while ( k >= 0 ) { j++; arg[j] = 1; k--; }
1340  j++;
1341  while ( j < i ) { arg[j] = 0; j++; }
1342  }
1343 #else
1344  j = i - 1;
1345  k = 0;
1346  while ( arg[j] == 1 && j >= 0 ) { j--; k++; }
1347  while ( arg[j] == 0 && j >= 0 ) j--;
1348  if ( j < 0 ) break;
1349  arg[j] = 0;
1350  while ( k >= 0 ) { j++; arg[j] = 1; k--; }
1351  j++;
1352  while ( j < i ) { arg[j] = 0; j++; }
1353 #endif
1354  }
1355  } while ( ntype == 0 && ++n <= i );
1356  AT.WorkPointer = oldwork;
1357  return(0);
1358 }
1359 
1360 /*
1361  #] DoDistrib :
1362  #[ EqualArg :
1363 
1364  Returns 1 if the arguments in the field are identical.
1365 */
1366 
1367 WORD EqualArg(WORD *parms, WORD num1, WORD num2)
1368 {
1369  WORD *t1, *t2;
1370  WORD i;
1371  t1 = parms;
1372  while ( --num1 >= 0 ) { NEXTARG(t1); }
1373  t2 = parms;
1374  while ( --num2 >= 0 ) { NEXTARG(t2); }
1375  if ( *t1 != *t2 ) return(0);
1376  if ( *t1 < 0 ) {
1377  if ( *t1 <= -FUNCTION || t1[1] == t2[1] ) return(1);
1378  return(0);
1379  }
1380  i = *t1;
1381  while ( --i >= 0 ) {
1382  if ( *t1 != *t2 ) return(0);
1383  t1++; t2++;
1384  }
1385  return(1);
1386 }
1387 
1388 /*
1389  #] EqualArg :
1390  #[ DoDelta3 :
1391 */
1392 
1393 WORD DoDelta3(PHEAD WORD *term, WORD level)
1394 {
1395  GETBIDENTITY
1396  WORD *t, *m, *m1, *m2, *stopper, *tstop, *termout, *dels, *taken;
1397  WORD *ic, *jc, *factors;
1398  WORD num, num2, i, j, k, knum, a;
1399  AN.TeInFun = AR.TePos = 0;
1400  tstop = term + *term;
1401  stopper = tstop - ABS(tstop[-1]);
1402  t = term+1;
1403  while ( ( *t != DELTA3 || ((t[1]-FUNHEAD) & 1 ) != 0 ) && t < stopper )
1404  t += t[1];
1405  if ( t >= stopper ) {
1406  MLOCK(ErrorMessageLock);
1407  MesPrint("Internal error with dd_ function");
1408  MUNLOCK(ErrorMessageLock);
1409  Terminate(-1);
1410  }
1411  m1 = t; m2 = t + t[1];
1412  num = t[1] - FUNHEAD;
1413  if ( num == 0 ) {
1414  termout = t = AT.WorkPointer;
1415  m = term;
1416  while ( m < m1 ) *t++ = *m++;
1417  m = m2; while ( m < tstop ) *t++ = *m++;
1418  *termout = WORDDIF(t,termout);
1419  AT.WorkPointer = t;
1420  if ( Generator(BHEAD termout,level) ) {
1421  MLOCK(ErrorMessageLock);
1422  MesCall("Do dd_");
1423  MUNLOCK(ErrorMessageLock);
1424  SETERROR(-1)
1425  }
1426  AT.WorkPointer = termout;
1427  return(0);
1428  }
1429  t += FUNHEAD;
1430 /*
1431  Step 1: sort the arguments
1432 */
1433  for ( i = 1; i < num; i++ ) {
1434  if ( t[i] < t[i-1] ) {
1435  a = t[i]; t[i] = t[i-1]; t[i-1] = a;
1436  j = i - 1;
1437  while ( j > 0 ) {
1438  if ( t[j] >= t[j-1] ) break;
1439  a = t[j]; t[j] = t[j-1]; t[j-1] = a;
1440  j--;
1441  }
1442  }
1443  }
1444 /*
1445  Step 2: Order them by occurrence
1446  In 'taken' we have the array with the number of occurrences.
1447  in 'dels' is the type of object.
1448 */
1449  m = taken = AT.WorkPointer;
1450  for ( i = 0; i < num; i++ ) *m++ = 0;
1451  dels = m; knum = 0;
1452  for ( i = 0; i < num; knum++ ) {
1453  *m++ = t[i]; i++; taken[knum] = 1;
1454  while ( i < num ) {
1455  if ( t[i] != t[i-1] ) break;
1456  i++; (taken[knum])++;
1457  }
1458  }
1459  for ( i = 0; i < knum; i++ ) *m++ = taken[i];
1460  ic = m; num2 = num/2;
1461  jc = ic + num2;
1462  factors = jc + num2;
1463  termout = factors + num2;
1464 /*
1465  The recursion has num/2 steps
1466 */
1467  k = 0;
1468  while ( k >= 0 ) {
1469  if ( k >= num2 ) {
1470  t = termout; m = term;
1471  while ( m < m1 ) *t++ = *m++;
1472  *t++ = DELTA; *t++ = num+2;
1473  for ( i = 0; i < num2; i++ ) {
1474  *t++ = dels[ic[i]]; *t++ = dels[jc[i]];
1475  }
1476  for ( i = 0; i < num2; i++ ) {
1477  if ( ic[i] == jc[i] ) {
1478  j = 1;
1479  while ( i < num2-1 && ic[i] == ic[i+1] && ic[i] == jc[i+1] )
1480  { i++; j++; }
1481  for ( a = 1; a < j; a++ ) {
1482  *t++ = SNUMBER; *t++ = 4; *t++ = 2*a+1; *t++ = 1;
1483  }
1484  for ( a = 0; a+1+i < num2; a++ ) {
1485  if ( ic[a+i] != ic[a+i+1] ) break;
1486  }
1487  if ( a > 0 ) {
1488  if ( GetBinom((UWORD *)(t+3),t+2,2*j+a,a) ) {
1489  MLOCK(ErrorMessageLock);
1490  MesCall("Do dd_");
1491  MUNLOCK(ErrorMessageLock);
1492  SETERROR(-1)
1493  }
1494  t[1] = ( t[2] < 0 ? -t[2]: t[2] ) + 3;
1495  *t = LNUMBER;
1496  t += t[1];
1497  }
1498  }
1499  else if ( factors[i] != 1 ) {
1500  *t++ = SNUMBER; *t++ = 4; *t++ = factors[i]; *t++ = 1;
1501  }
1502  }
1503  for ( i = 0; i < num2-1; i++ ) {
1504  if ( ic[i] == jc[i] ) continue;
1505  j = 1;
1506  while ( i < num2-1 && jc[i] == jc[i+1] && ic[i] == ic[i+1] ) {
1507  i++; j++;
1508  }
1509  for ( a = 0; a+i < num2-1; a++ ) {
1510  if ( ic[i+a] != ic[i+a+1] ) break;
1511  }
1512  if ( a > 0 ) {
1513  if ( GetBinom((UWORD *)(t+3),t+2,j+a,a) ) {
1514  MLOCK(ErrorMessageLock);
1515  MesCall("Do dd_");
1516  MUNLOCK(ErrorMessageLock);
1517  SETERROR(-1)
1518  }
1519  t[1] = ( t[2] < 0 ? -t[2]: t[2] ) + 3;
1520  *t = LNUMBER;
1521  t += t[1];
1522  }
1523  }
1524  m = m2;
1525  while ( m < tstop ) *t++ = *m++;
1526  *termout = WORDDIF(t,termout);
1527  AT.WorkPointer = t;
1528  if ( Generator(BHEAD termout,level) ) {
1529  MLOCK(ErrorMessageLock);
1530  MesCall("Do dd_");
1531  MUNLOCK(ErrorMessageLock);
1532  SETERROR(-1)
1533  }
1534  k--;
1535  if ( k >= 0 ) goto nextj;
1536  else break;
1537  }
1538  for ( ic[k] = 0; ic[k] < knum; ic[k]++ ) {
1539  if ( taken[ic[k]] > 0 ) break;
1540  }
1541  if ( k > 0 && ic[k-1] == ic[k] ) jc[k] = jc[k-1];
1542  else jc[k] = ic[k];
1543  for ( ; jc[k] < knum; jc[k]++ ) {
1544  if ( taken[jc[k]] <= 0 ) continue;
1545  if ( ic[k] == jc[k] ) {
1546  if ( taken[jc[k]] <= 1 ) continue;
1547 /*
1548  factors[k] = taken[ic[k]];
1549  if ( ( factors[k] & 1 ) == 0 ) (factors[k])--;
1550 */
1551  taken[ic[k]] -= 2;
1552  }
1553  else {
1554  factors[k] = taken[jc[k]];
1555  (taken[ic[k]])--; (taken[jc[k]])--;
1556  }
1557  k++;
1558  goto nextk; /* This is the simulated recursion */
1559 nextj:;
1560  (taken[ic[k]])++; (taken[jc[k]])++;
1561  }
1562  k--;
1563  if ( k >= 0 ) goto nextj;
1564 nextk:;
1565  }
1566  AT.WorkPointer = taken;
1567  return(0);
1568 }
1569 
1570 /*
1571  #] DoDelta3 :
1572  #] Distribute :
1573  #[ DoShuffle :
1574 
1575  Merges the arguments of all occurrences of function fun into a
1576  single occurrence of fun. The opposite of Distrib_
1577  Syntax:
1578  Shuffle[,once|all],fun;
1579  Shuffle[,once|all],$fun;
1580  The expansion of the dollar should give a single function.
1581  The dollar is indicated as usual with a negative value.
1582  option = 1 (once): generate identical results only once
1583  option = 0 (all): generate identical results with combinatorics (default)
1584 */
1585 
1586 /*
1587  We use the Shuffle routine which has a large amount of combinatorics.
1588  It doesn't have grouped combinatorics as in (0,1,2)*(0,1,3) where the
1589  groups (0,1) also cause double terms.
1590 */
1591 
1592 WORD DoShuffle(PHEAD WORD *term, WORD level, WORD fun, WORD option)
1593 {
1594  GETBIDENTITY
1595  SHvariables SHback, *SH = &(AN.SHvar);
1596  WORD *t1, *t2, *tstop, ncoef, n = fun, *to, *from;
1597  int i, error;
1598  LONG k;
1599  UWORD *newcombi;
1600 
1601  if ( n < 0 ) {
1602  if ( ( n = DolToFunction(BHEAD -n) ) == 0 ) {
1603  MLOCK(ErrorMessageLock);
1604  MesPrint("$-variable in merge statement did not evaluate to a function.");
1605  MUNLOCK(ErrorMessageLock);
1606  return(1);
1607  }
1608  }
1609  if ( AT.WorkPointer + 3*(*term) + AM.MaxTal > AT.WorkTop ) {
1610  MLOCK(ErrorMessageLock);
1611  MesWork();
1612  MUNLOCK(ErrorMessageLock);
1613  return(-1);
1614  }
1615 
1616  tstop = term + *term;
1617  ncoef = tstop[-1];
1618  tstop -= ABS(ncoef);
1619  t1 = term + 1;
1620  while ( t1 < tstop ) {
1621  if ( ( *t1 == n ) && ( t1+t1[1] < tstop ) && ( t1[1] > FUNHEAD ) ) {
1622  t2 = t1 + t1[1];
1623  if ( t2 >= tstop ) {
1624  return(Generator(BHEAD term,level));
1625  }
1626  while ( t2 < tstop ) {
1627  if ( ( *t2 == n ) && ( t2[1] > FUNHEAD ) ) break;
1628  t2 += t2[1];
1629  }
1630  if ( t2 < tstop ) break;
1631  }
1632  t1 += t1[1];
1633  }
1634  if ( t1 >= tstop ) {
1635  return(Generator(BHEAD term,level));
1636  }
1637  *AN.RepPoint = 1;
1638 /*
1639  Now we have two occurrences of the function.
1640  Back up all relevant variables and load all the stuff that needs to be
1641  passed on.
1642 */
1643  SHback = AN.SHvar;
1644  SH->finishuf = &FinishShuffle;
1645  SH->do_uffle = &DoShuffle;
1646  SH->outterm = AT.WorkPointer;
1647  AT.WorkPointer += *term;
1648  SH->stop1 = t1 + t1[1];
1649  SH->stop2 = t2 + t2[1];
1650  SH->thefunction = n;
1651  SH->option = option;
1652  SH->level = level;
1653  SH->incoef = tstop;
1654  SH->nincoef = ncoef;
1655 
1656  if ( AN.SHcombi == 0 || AN.SHcombisize == 0 ) {
1657  AN.SHcombisize = 200;
1658  AN.SHcombi = (UWORD *)Malloc1(AN.SHcombisize*sizeof(UWORD),"AN.SHcombi");
1659  SH->combilast = 0;
1660  SHback.combilast = 0;
1661  }
1662  else {
1663  SH->combilast += AN.SHcombi[SH->combilast]+1;
1664  if ( SH->combilast >= AN.SHcombisize - 100 ) {
1665  newcombi = (UWORD *)Malloc1(2*AN.SHcombisize*sizeof(UWORD),"AN.SHcombi");
1666  for ( k = 0; k < AN.SHcombisize; k++ ) newcombi[k] = AN.SHcombi[k];
1667  M_free(AN.SHcombi,"AN.SHcombi");
1668  AN.SHcombi = newcombi;
1669  AN.SHcombisize *= 2;
1670  }
1671  }
1672  AN.SHcombi[SH->combilast] = 1;
1673  AN.SHcombi[SH->combilast+1] = 1;
1674 
1675  i = t1-term; to = SH->outterm; from = term;
1676  NCOPY(to,from,i)
1677  SH->outfun = to;
1678  for ( i = 0; i < FUNHEAD; i++ ) { *to++ = t1[i]; }
1679 
1680  error = Shuffle(BHEAD t1+FUNHEAD,t2+FUNHEAD,to);
1681 
1682  AT.WorkPointer = SH->outterm;
1683  AN.SHvar = SHback;
1684  if ( error ) {
1685  MesCall("DoShuffle");
1686  return(-1);
1687  }
1688  return(0);
1689 }
1690 
1691 /*
1692  #] DoShuffle :
1693  #[ Shuffle :
1694 
1695  How to make shuffles:
1696 
1697  We have two lists of arguments. We have to make a single
1698  shuffle of them. All combinations. Doubles should have as
1699  much as possible a combinatorics factor. Sometimes this is
1700  very difficult as in:
1701  (0,1,2)x(0,1,3) = -> (0,1) is a repeated pattern and the
1702  factor on that is difficult
1703  Simple way: (without combinatorics)
1704  repeat id f0(?c)*f(x1?,?a)*f(x2?,?b) =
1705  +f0(?c,x1)*f(?a)*f(x2,?b)
1706  +f0(?c,x2)*f(x1,?a)*f(?b);
1707  Refinement:
1708  if ( x1 == x2 ) check how many more there are of the same.
1709  --> (n1,x) and (n2,x)
1710  id f0(?c)*f1((n1,x),?b)*f2((n2,x),?c) =
1711  +binom_(n1+n2,n1)*f0(?c,(n1+n2,x))*f1(?a)*f2(?b)
1712  +sum_(j,0,n1-1,binom_(n2+j,j)*f0(?c,(j+n2,x))
1713  *f1((n1-j),?a)*f2(?b))*force2
1714  +sum_(j,0,n2-1,binom_(n1+j,j)*f0(?c,(j+n1,x))
1715  *f1(?a)*f2((n2-j),?b))*force1
1716  The force operation can be executed directly
1717 
1718  The next question is how to program this: recursively or linearly
1719  which would require simulation of a recursion. Recursive is clearest
1720  but we need to pass a number of arguments from the calling routine
1721  to the final routine. This is done with AN.SHvar.
1722 
1723  We need space for the accumulation of the combinatoric factors.
1724 */
1725 
1726 int Shuffle(PHEAD WORD *from1, WORD *from2, WORD *to)
1727 {
1728  WORD *t, *fr, *next1, *next2, na, *fn1, *fn2, *tt;
1729  int i, n, n1, n2, j;
1730  LONG combilast;
1731  SHvariables *SH = &(AN.SHvar);
1732  if ( from1 == SH->stop1 && from2 == SH->stop2 ) {
1733  return(FiniShuffle(BHEAD to));
1734  }
1735  else if ( from1 == SH->stop1 ) {
1736  i = SH->stop2 - from2; t = to; tt = from2; NCOPY(t,tt,i)
1737  return(FiniShuffle(BHEAD t));
1738  }
1739  else if ( from2 == SH->stop2 ) {
1740  i = SH->stop1 - from1; t = to; tt = from1; NCOPY(t,tt,i)
1741  return(FiniShuffle(BHEAD t));
1742  }
1743 /*
1744  Compare lead arguments
1745 */
1746  if ( AreArgsEqual(from1,from2) ) {
1747 /*
1748  First find out how many of each
1749 */
1750  next1 = from1; n1 = 1; NEXTARG(next1)
1751  while ( ( next1 < SH->stop1 ) && AreArgsEqual(from1,next1) ) {
1752  n1++; NEXTARG(next1)
1753  }
1754  next2 = from2; n2 = 1; NEXTARG(next2)
1755  while ( ( next2 < SH->stop2 ) && AreArgsEqual(from2,next2) ) {
1756  n2++; NEXTARG(next2)
1757  }
1758  combilast = SH->combilast;
1759 /*
1760  +binom_(n1+n2,n1)*f0(?c,(n1+n2,x))*f1(?a)*f2(?b)
1761 */
1762  t = to;
1763  n = n1 + n2;
1764  while ( --n >= 0 ) { fr = from1; CopyArg(t,fr) }
1765  if ( GetBinom((UWORD *)(t),&na,n1+n2,n1) ) goto shuffcall;
1766  if ( combilast + AN.SHcombi[combilast] + na + 2 >= AN.SHcombisize ) {
1767 /*
1768  We need more memory in this stack. Fortunately this is the
1769  only place where we have to do this, because the other factors
1770  are definitely smaller.
1771  Layout: size, LongInteger, size, LongInteger, .....
1772  We start pointing at the last one.
1773 */
1774  UWORD *combi = (UWORD *)Malloc1(2*AN.SHcombisize*2,"AN.SHcombi");
1775  LONG jj;
1776  for ( jj = 0; jj < AN.SHcombisize; jj++ ) combi[jj] = AN.SHcombi[jj];
1777  AN.SHcombisize *= 2;
1778  M_free(AN.SHcombi,"AN.SHcombi");
1779  AN.SHcombi = combi;
1780  }
1781  if ( MulLong((UWORD *)(AN.SHcombi+combilast+1),AN.SHcombi[combilast],
1782  (UWORD *)(t),na,
1783  (UWORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+2),
1784  (WORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+1)) ) goto shuffcall;
1785  SH->combilast = combilast + AN.SHcombi[combilast] + 1;
1786  if ( next1 >= SH->stop1 ) {
1787  fr = next2; i = SH->stop2 - fr;
1788  NCOPY(t,fr,i)
1789  if ( FiniShuffle(BHEAD t) ) goto shuffcall;
1790  }
1791  else if ( next2 >= SH->stop2 ) {
1792  fr = next1; i = SH->stop1 - fr;
1793  NCOPY(t,fr,i)
1794  if ( FiniShuffle(BHEAD t) ) goto shuffcall;
1795  }
1796  else {
1797  if ( Shuffle(BHEAD next1,next2,t) ) goto shuffcall;
1798  }
1799  SH->combilast = combilast;
1800 /*
1801  +sum_(j,0,n1-1,binom_(n2+j,j)*f0(?c,(j+n2,x))
1802  *f1((n1-j),?a)*f2(?b))*force2
1803 */
1804  if ( next2 < SH->stop2 ) {
1805  t = to;
1806  n = n2;
1807  while ( --n >= 0 ) { fr = from1; CopyArg(t,fr) }
1808  for ( j = 0; j < n1; j++ ) {
1809  if ( GetBinom((UWORD *)(t),&na,n2+j,j) ) goto shuffcall;
1810  if ( MulLong((UWORD *)(AN.SHcombi+combilast+1),AN.SHcombi[combilast],
1811  (UWORD *)(t),na,
1812  (UWORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+2),
1813  (WORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+1)) ) goto shuffcall;
1814  SH->combilast = combilast + AN.SHcombi[combilast] + 1;
1815  if ( j > 0 ) { fr = from1; CopyArg(t,fr) }
1816  fn2 = next2; tt = t;
1817  CopyArg(tt,fn2)
1818 
1819  if ( fn2 >= SH->stop2 ) {
1820  n = n1-j;
1821  while ( --n >= 0 ) { fr = from1; CopyArg(tt,fr) }
1822  fr = next1; i = SH->stop1 - fr;
1823  NCOPY(tt,fr,i)
1824  if ( FiniShuffle(BHEAD tt) ) goto shuffcall;
1825  }
1826  else {
1827  n = j; fn1 = from1; while ( --n >= 0 ) { NEXTARG(fn1) }
1828  if ( Shuffle(BHEAD fn1,fn2,tt) ) goto shuffcall;
1829  }
1830  SH->combilast = combilast;
1831  }
1832  }
1833 /*
1834  +sum_(j,0,n2-1,binom_(n1+j,j)*f0(?c,(j+n1,x))
1835  *f1(?a)*f2((n2-j),?b))*force1
1836 */
1837  if ( next1 < SH->stop1 ) {
1838  t = to;
1839  n = n1;
1840  while ( --n >= 0 ) { fr = from1; CopyArg(t,fr) }
1841  for ( j = 0; j < n2; j++ ) {
1842  if ( GetBinom((UWORD *)(t),&na,n1+j,j) ) goto shuffcall;
1843  if ( MulLong((UWORD *)(AN.SHcombi+combilast+1),AN.SHcombi[combilast],
1844  (UWORD *)(t),na,
1845  (UWORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+2),
1846  (WORD *)(AN.SHcombi+combilast+AN.SHcombi[combilast]+1)) ) goto shuffcall;
1847  SH->combilast = combilast + AN.SHcombi[combilast] + 1;
1848  if ( j > 0 ) { fr = from1; CopyArg(t,fr) }
1849  fn1 = next1; tt = t;
1850  CopyArg(tt,fn1)
1851 
1852  if ( fn1 >= SH->stop1 ) {
1853  n = n2-j;
1854  while ( --n >= 0 ) { fr = from1; CopyArg(tt,fr) }
1855  fr = next2; i = SH->stop2 - fr;
1856  NCOPY(tt,fr,i)
1857  if ( FiniShuffle(BHEAD tt) ) goto shuffcall;
1858  }
1859  else {
1860  n = j; fn2 = from2; while ( --n >= 0 ) { NEXTARG(fn2) }
1861  if ( Shuffle(BHEAD fn1,fn2,tt) ) goto shuffcall;
1862  }
1863  SH->combilast = combilast;
1864  }
1865  }
1866  }
1867  else {
1868 /*
1869  Argument from first list
1870 */
1871  t = to;
1872  fr = from1;
1873  CopyArg(t,fr)
1874  if ( fr >= SH->stop1 ) {
1875  fr = from2; i = SH->stop2 - fr;
1876  NCOPY(t,fr,i)
1877  if ( FiniShuffle(BHEAD t) ) goto shuffcall;
1878  }
1879  else {
1880  if ( Shuffle(BHEAD fr,from2,t) ) goto shuffcall;
1881  }
1882 /*
1883  Argument from second list
1884 */
1885  t = to;
1886  fr = from2;
1887  CopyArg(t,fr)
1888  if ( fr >= SH->stop2 ) {
1889  fr = from1; i = SH->stop1 - fr;
1890  NCOPY(t,fr,i)
1891  if ( FiniShuffle(BHEAD t) ) goto shuffcall;
1892  }
1893  else {
1894  if ( Shuffle(BHEAD from1,fr,t) ) goto shuffcall;
1895  }
1896  }
1897  return(0);
1898 shuffcall:
1899  MesCall("Shuffle");
1900  return(-1);
1901 }
1902 
1903 /*
1904  #] Shuffle :
1905  #[ FinishShuffle :
1906 
1907  The complications here are:
1908  1: We want to save space. We put the output term in 'out' straight
1909  on top of what we produced thusfar. We have to copy the early
1910  piece because once the term goes back to Generator, Normalize can
1911  change it in situ
1912  2: There can be other occurrence of the function between the two
1913  that we did. For shuffles that isn't likely, but we use this
1914  routine also for the stuffles and there it can happen.
1915 */
1916 
1917 int FinishShuffle(PHEAD WORD *fini)
1918 {
1919  WORD *t, *t1, *oldworkpointer = AT.WorkPointer, *tcoef, ntcoef, *out;
1920  int i;
1921  SHvariables *SH = &(AN.SHvar);
1922  SH->outfun[1] = fini - SH->outfun;
1923  if ( functions[SH->outfun[0]-FUNCTION].symmetric != 0 )
1924  SH->outfun[2] |= DIRTYSYMFLAG;
1925  out = fini; i = fini - SH->outterm; t = SH->outterm;
1926  NCOPY(fini,t,i)
1927  t = SH->stop1;
1928  t1 = t + t[1];
1929  while ( t1 < SH->stop2 ) { t = t1; t1 = t + t[1]; }
1930  t1 = SH->stop1;
1931  while ( t1 < t ) *fini++ = *t1++;
1932  t = SH->stop2;
1933  while ( t < SH->incoef ) *fini++ = *t++;
1934  tcoef = fini;
1935  ntcoef = SH->nincoef;
1936  i = ABS(ntcoef);
1937  NCOPY(fini,t,i);
1938  ntcoef = REDLENG(ntcoef);
1939  Mully(BHEAD (UWORD *)tcoef,&ntcoef,
1940  (UWORD *)(AN.SHcombi+SH->combilast+1),AN.SHcombi[SH->combilast]);
1941  ntcoef = INCLENG(ntcoef);
1942  fini = tcoef + ABS(ntcoef);
1943  if ( ( ( SH->option & 2 ) != 0 ) && ( ( SH->option & 256 ) != 0 ) ) ntcoef = -ntcoef;
1944  fini[-1] = ntcoef;
1945  i = *out = fini - out;
1946 /*
1947  Now check whether we have to do more
1948 */
1949  AT.WorkPointer = out + *out;
1950  if ( ( SH->option & 1 ) == 1 ) {
1951  if ( Generator(BHEAD out,SH->level) ) goto Finicall;
1952  }
1953  else {
1954  if ( DoShtuffle(BHEAD out,SH->level,SH->thefunction,SH->option) ) goto Finicall;
1955  }
1956  AT.WorkPointer = oldworkpointer;
1957  return(0);
1958 Finicall:
1959  AT.WorkPointer = oldworkpointer;
1960  MesCall("FinishShuffle");
1961  return(-1);
1962 }
1963 
1964 /*
1965  #] FinishShuffle :
1966  #[ DoStuffle :
1967 
1968  Stuffling is a variation of shuffling.
1969  In the stuffling we insist that the arguments are (short) integers. nonzero.
1970  The stuffle sum is x st y = sig_(x)*sig_(y)*(abs(x)+abs(y))
1971  The way we do this is:
1972  1: count the arguments in each function: n1, n2
1973  2: take the minimum minval = min(n1,n2).
1974  3: for ( j = 0; j <= min; j++ ) take j elements in each of the lists.
1975  4: the j+1 groups of remaining arguments have to each be shuffled
1976  5: the j selected pairs have to be stuffle added.
1977  We can use many of the shuffle things.
1978  Considering the recursive nature of the generation we actually don't
1979  need to know n1, n2, minval.
1980 */
1981 
1982 WORD DoStuffle(PHEAD WORD *term, WORD level, WORD fun, WORD option)
1983 {
1984  GETBIDENTITY
1985  SHvariables SHback, *SH = &(AN.SHvar);
1986  WORD *t1, *t2, *tstop, *t1stop, *t2stop, ncoef, n = fun, *to, *from;
1987  WORD *r1, *r2;
1988  int i, error;
1989  LONG k;
1990  UWORD *newcombi;
1991 #ifdef NEWCODE
1992  WORD *rr1, *rr2, i1, i2;
1993 #endif
1994  if ( n < 0 ) {
1995  if ( ( n = DolToFunction(BHEAD -n) ) == 0 ) {
1996  MLOCK(ErrorMessageLock);
1997  MesPrint("$-variable in merge statement did not evaluate to a function.");
1998  MUNLOCK(ErrorMessageLock);
1999  return(1);
2000  }
2001  }
2002  if ( AT.WorkPointer + 3*(*term) + AM.MaxTal > AT.WorkTop ) {
2003  MLOCK(ErrorMessageLock);
2004  MesWork();
2005  MUNLOCK(ErrorMessageLock);
2006  return(-1);
2007  }
2008 
2009  tstop = term + *term;
2010  ncoef = tstop[-1];
2011  tstop -= ABS(ncoef);
2012  t1 = term + 1;
2013 retry1:;
2014  while ( t1 < tstop ) {
2015  if ( ( *t1 == n ) && ( t1+t1[1] < tstop ) && ( t1[1] > FUNHEAD ) ) {
2016  t2 = t1 + t1[1];
2017  if ( t2 >= tstop ) {
2018  return(Generator(BHEAD term,level));
2019  }
2020 retry2:;
2021  while ( t2 < tstop ) {
2022  if ( ( *t2 == n ) && ( t2[1] > FUNHEAD ) ) break;
2023  t2 += t2[1];
2024  }
2025  if ( t2 < tstop ) break;
2026  }
2027  t1 += t1[1];
2028  }
2029  if ( t1 >= tstop ) {
2030  return(Generator(BHEAD term,level));
2031  }
2032 /*
2033  Next we have to check that the arguments are of the correct type
2034  At the same time we can count them.
2035 */
2036 #ifndef NEWCODE
2037  t1stop = t1 + t1[1];
2038  r1 = t1 + FUNHEAD;
2039  while ( r1 < t1stop ) {
2040  if ( *r1 != -SNUMBER ) break;
2041  if ( r1[1] == 0 ) break;
2042  r1 += 2;
2043  }
2044  if ( r1 < t1stop ) { t1 = t2; goto retry1; }
2045  t2stop = t2 + t2[1];
2046  r2 = t2 + FUNHEAD;
2047  while ( r2 < t2stop ) {
2048  if ( *r2 != -SNUMBER ) break;
2049  if ( r2[1] == 0 ) break;
2050  r2 += 2;
2051  }
2052  if ( r2 < t2stop ) { t2 = t2 + t2[1]; goto retry2; }
2053 #else
2054  t1stop = t1 + t1[1];
2055  r1 = t1 + FUNHEAD;
2056  while ( r1 < t1stop ) {
2057  if ( *r1 == -SNUMBER ) {
2058  if ( r1[1] == 0 ) break;
2059  r1 += 2; continue;
2060  }
2061  else if ( *r1 == -SYMBOL ) {
2062  if ( ( symbols[r1[1]].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
2063  break;
2064  r1 += 2; continue;
2065  }
2066  if ( *r1 > 0 && *r1 == r1[ARGHEAD]+ARGHEAD ) {
2067  if ( ABS(r1[r1[0]-1]) == r1[0]-ARGHEAD-1 ) {}
2068  else if ( r1[ARGHEAD+1] == SYMBOL ) {
2069  rr1 = r1 + ARGHEAD + 3;
2070  i1 = rr1[-1]-2;
2071  while ( i1 > 0 ) {
2072  if ( ( symbols[*rr1].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
2073  break;
2074  i1 -= 2; rr1 += 2;
2075  }
2076  if ( i1 > 0 ) break;
2077  }
2078  else break;
2079  rr1 = r1+*r1-1;
2080  i1 = (ABS(*rr1)-1)/2;
2081  while ( i1 > 1 ) {
2082  if ( rr1[-1] ) break;
2083  i1--; rr1--;
2084  }
2085  if ( i1 > 1 || rr1[-1] != 1 ) break;
2086  r1 += *r1;
2087  }
2088  else break;
2089  }
2090  if ( r1 < t1stop ) { t1 = t2; goto retry1; }
2091  t2stop = t2 + t2[1];
2092  r2 = t2 + FUNHEAD;
2093 
2094  while ( r2 < t2stop ) {
2095  if ( *r2 == -SNUMBER ) {
2096  if ( r2[1] == 0 ) break;
2097  r2 += 2; continue;
2098  }
2099  else if ( *r2 == -SYMBOL ) {
2100  if ( ( symbols[r2[1]].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
2101  break;
2102  r2 += 2; continue;
2103  }
2104  if ( *r2 > 0 && *r2 == r2[ARGHEAD]+ARGHEAD ) {
2105  if ( ABS(r2[r2[0]-1]) == r2[0]-ARGHEAD-1 ) {}
2106  else if ( r2[ARGHEAD+1] == SYMBOL ) {
2107  rr2 = r2 + ARGHEAD + 3;
2108  i2 = rr2[-1]-2;
2109  while ( i2 > 0 ) {
2110  if ( ( symbols[*rr2].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
2111  break;
2112  i2 -= 2; rr2 += 2;
2113  }
2114  if ( i2 > 0 ) break;
2115  }
2116  else break;
2117  rr2 = r2+*r2-1;
2118  i2 = (ABS(*rr2)-1)/2;
2119  while ( i2 > 1 ) {
2120  if ( rr2[-1] ) break;
2121  i2--; rr2--;
2122  }
2123  if ( i2 > 1 || rr2[-1] != 1 ) break;
2124  r2 += *r2;
2125  }
2126  else break;
2127  }
2128  if ( r2 < t2stop ) { t2 = t2 + t2[1]; goto retry2; }
2129 #endif
2130 /*
2131  OK, now we got two objects that can be used.
2132 */
2133  *AN.RepPoint = 1;
2134 
2135  SHback = AN.SHvar;
2136  SH->finishuf = &FinishStuffle;
2137  SH->do_uffle = &DoStuffle;
2138  SH->outterm = AT.WorkPointer;
2139  AT.WorkPointer += *term;
2140  SH->ststop1 = t1 + t1[1];
2141  SH->ststop2 = t2 + t2[1];
2142  SH->thefunction = n;
2143  SH->option = option;
2144  SH->level = level;
2145  SH->incoef = tstop;
2146  SH->nincoef = ncoef;
2147  if ( AN.SHcombi == 0 || AN.SHcombisize == 0 ) {
2148  AN.SHcombisize = 200;
2149  AN.SHcombi = (UWORD *)Malloc1(AN.SHcombisize*sizeof(UWORD),"AN.SHcombi");
2150  SH->combilast = 0;
2151  SHback.combilast = 0;
2152  }
2153  else {
2154  SH->combilast += AN.SHcombi[SH->combilast]+1;
2155  if ( SH->combilast >= AN.SHcombisize - 100 ) {
2156  newcombi = (UWORD *)Malloc1(2*AN.SHcombisize*sizeof(UWORD),"AN.SHcombi");
2157  for ( k = 0; k < AN.SHcombisize; k++ ) newcombi[k] = AN.SHcombi[k];
2158  M_free(AN.SHcombi,"AN.SHcombi");
2159  AN.SHcombi = newcombi;
2160  AN.SHcombisize *= 2;
2161  }
2162  }
2163  AN.SHcombi[SH->combilast] = 1;
2164  AN.SHcombi[SH->combilast+1] = 1;
2165 
2166  i = t1-term; to = SH->outterm; from = term;
2167  NCOPY(to,from,i)
2168  SH->outfun = to;
2169  for ( i = 0; i < FUNHEAD; i++ ) { *to++ = t1[i]; }
2170 
2171  error = Stuffle(BHEAD t1+FUNHEAD,t2+FUNHEAD,to);
2172 
2173  AT.WorkPointer = SH->outterm;
2174  AN.SHvar = SHback;
2175  if ( error ) {
2176  MesCall("DoStuffle");
2177  return(-1);
2178  }
2179  return(0);
2180 }
2181 
2182 /*
2183  #] DoStuffle :
2184  #[ Stuffle :
2185 
2186  The way to generate the stuffles
2187  1: select an argument in the first list (for(j1=0;j1<last;j1++))
2188  2: select an argument in the second list (for(j2=0;j2<last;j2++))
2189  3: put values for SH->ststop1 and SH->ststop2 at these arguments.
2190  4: generate all shuffles of the arguments in front.
2191  5: Then put the stuffle sum of arg(j1) and arg(j2)
2192  6: Then continue calling Stuffle
2193  7: Once one gets exhausted, we can clean up the list and call FinishShuffle
2194  8: if ( ( SH->option & 2 ) != 0 ) the stuffle sum is negative.
2195 */
2196 
2197 int Stuffle(PHEAD WORD *from1, WORD *from2, WORD *to)
2198 {
2199  GETBIDENTITY
2200  WORD *t, *tf, *next1, *next2, *st1, *st2, *save1, *save2;
2201  SHvariables *SH = &(AN.SHvar);
2202  int i, retval;
2203 /*
2204  First the special cases (exhausted list(s)):
2205 */
2206  save1 = SH->stop1; save2 = SH->stop2;
2207  if ( from1 >= SH->ststop1 && from2 == SH->ststop2 ) {
2208  SH->stop1 = SH->ststop1;
2209  SH->stop2 = SH->ststop2;
2210  retval = FinishShuffle(BHEAD to);
2211  SH->stop1 = save1; SH->stop2 = save2;
2212  return(retval);
2213  }
2214  else if ( from1 >= SH->ststop1 ) {
2215  i = SH->ststop2 - from2; t = to; tf = from2; NCOPY(t,tf,i)
2216  SH->stop1 = SH->ststop1;
2217  SH->stop2 = SH->ststop2;
2218  retval = FinishShuffle(BHEAD t);
2219  SH->stop1 = save1; SH->stop2 = save2;
2220  return(retval);
2221  }
2222  else if ( from2 >= SH->ststop2 ) {
2223  i = SH->ststop1 - from1; t = to; tf = from1; NCOPY(t,tf,i)
2224  SH->stop1 = SH->ststop1;
2225  SH->stop2 = SH->ststop2;
2226  retval = FinishShuffle(BHEAD t);
2227  SH->stop1 = save1; SH->stop2 = save2;
2228  return(retval);
2229  }
2230 /*
2231  Now the case that we have no stuffle sums.
2232 */
2233  SH->stop1 = SH->ststop1;
2234  SH->stop2 = SH->ststop2;
2235  SH->finishuf = &FinishShuffle;
2236  if ( Shuffle(BHEAD from1,from2,to) ) goto stuffcall;
2237  SH->finishuf = &FinishStuffle;
2238 /*
2239  Now we have to select a pair, one from 1 and one from 2.
2240 */
2241 #ifndef NEWCODE
2242  st1 = from1; next1 = st1+2; /* <----- */
2243 #else
2244  st1 = next1 = from1;
2245  NEXTARG(next1)
2246 #endif
2247  while ( next1 <= SH->ststop1 ) {
2248 #ifndef NEWCODE
2249  st2 = from2; next2 = st2+2; /* <----- */
2250 #else
2251  next2 = st2 = from2;
2252  NEXTARG(next2)
2253 #endif
2254  while ( next2 <= SH->ststop2 ) {
2255  SH->stop1 = st1;
2256  SH->stop2 = st2;
2257  if ( st1 == from1 && st2 == from2 ) {
2258  t = to;
2259 #ifndef NEWCODE
2260  *t++ = -SNUMBER; *t++ = StuffAdd(st1[1],st2[1]);
2261 #else
2262  t = StuffRootAdd(st1,st2,t);
2263 #endif
2264  SH->option ^= 256;
2265  if ( Stuffle(BHEAD next1,next2,t) ) goto stuffcall;
2266  SH->option ^= 256;
2267  }
2268  else if ( st1 == from1 ) {
2269  i = st2-from2;
2270  t = to; tf = from2; NCOPY(t,tf,i)
2271 #ifndef NEWCODE
2272  *t++ = -SNUMBER; *t++ = StuffAdd(st1[1],st2[1]);
2273 #else
2274  t = StuffRootAdd(st1,st2,t);
2275 #endif
2276  SH->option ^= 256;
2277  if ( Stuffle(BHEAD next1,next2,t) ) goto stuffcall;
2278  SH->option ^= 256;
2279  }
2280  else if ( st2 == from2 ) {
2281  i = st1-from1;
2282  t = to; tf = from1; NCOPY(t,tf,i)
2283 #ifndef NEWCODE
2284  *t++ = -SNUMBER; *t++ = StuffAdd(st1[1],st2[1]);
2285 #else
2286  t = StuffRootAdd(st1,st2,t);
2287 #endif
2288  SH->option ^= 256;
2289  if ( Stuffle(BHEAD next1,next2,t) ) goto stuffcall;
2290  SH->option ^= 256;
2291  }
2292  else {
2293  if ( Shuffle(BHEAD from1,from2,to) ) goto stuffcall;
2294  }
2295 #ifndef NEWCODE
2296  st2 = next2; next2 += 2; /* <----- */
2297 #else
2298  st2 = next2;
2299  NEXTARG(next2)
2300 #endif
2301  }
2302 #ifndef NEWCODE
2303  st1 = next1; next1 += 2; /* <----- */
2304 #else
2305  st1 = next1;
2306  NEXTARG(next1)
2307 #endif
2308  }
2309  SH->stop1 = save1; SH->stop2 = save2;
2310  return(0);
2311 stuffcall:;
2312  MesCall("Stuffle");
2313  return(-1);
2314 }
2315 
2316 /*
2317  #] Stuffle :
2318  #[ FinishStuffle :
2319 
2320  The program only comes here from the Shuffle routine.
2321  It should add the stuffle sum and then call Stuffle again.
2322 */
2323 
2324 int FinishStuffle(PHEAD WORD *fini)
2325 {
2326  GETBIDENTITY
2327  SHvariables *SH = &(AN.SHvar);
2328 #ifdef NEWCODE
2329  WORD *next1 = SH->stop1, *next2 = SH->stop2;
2330  fini = StuffRootAdd(next1,next2,fini);
2331 #else
2332  *fini++ = -SNUMBER; *fini++ = StuffAdd(SH->stop1[1],SH->stop2[1]);
2333 #endif
2334  SH->option ^= 256;
2335 #ifdef NEWCODE
2336  NEXTARG(next1)
2337  NEXTARG(next2)
2338  if ( Stuffle(BHEAD next1,next2,fini) ) goto stuffcall;
2339 #else
2340  if ( Stuffle(BHEAD SH->stop1+2,SH->stop2+2,fini) ) goto stuffcall;
2341 #endif
2342  SH->option ^= 256;
2343  return(0);
2344 stuffcall:;
2345  MesCall("FinishStuffle");
2346  return(-1);
2347 }
2348 
2349 /*
2350  #] FinishStuffle :
2351  #[ StuffRootAdd :
2352 
2353  Makes the stuffle sum of two arguments.
2354  The arguments can be of one of three types:
2355  1: -SNUMBER,num
2356  2: -SYMBOL,symbol
2357  3: Numerical (long) argument.
2358  4: Generic argument with (only) symbols that are roots of unity and
2359  a coefficient.
2360  We have excluded the case that both t1 and t2 are of type 1:
2361  The output should be written to 'to' and the new fill position should
2362  be the return value.
2363  `to' is inside the workspace.
2364 
2365  The stuffle sum is sig_(t2)*t1+sig_(t1)*t2
2366  or sig_(t1)*sig_(t2)*(abs_(t1)+abs_(t2))
2367 */
2368 
2369 #ifdef NEWCODE
2370 
2371 WORD *StuffRootAdd(WORD *t1, WORD *t2, WORD *to)
2372 {
2373  int type1, type2, type3, sgn, sgn1, sgn2, sgn3, pow, root, nosymbols, i;
2374  WORD *tt1, *tt2, it1, it2, *t3, *r, size1, size2, size3;
2375  WORD scratch[2];
2376  LONG x;
2377  if ( *t1 == -SNUMBER ) { type1 = 1; if ( t1[1] < 0 ) sgn1 = -1; else sgn1 = 1; }
2378  else if ( *t1 == -SYMBOL ) { type1 = 2; sgn1 = 1; }
2379  else if ( ABS(t1[*t1-1]) == *t1-ARGHEAD-1 ) {
2380  type1 = 3; if ( t1[*t1-1] < 0 ) sgn1 = -1; else sgn1 = 1; }
2381  else { type1 = 4; if ( t1[*t1-1] < 0 ) sgn1 = -1; else sgn1 = 1; }
2382  if ( *t2 == -SNUMBER ) { type2 = 1; if ( t2[1] < 0 ) sgn2 = -1; else sgn2 = 1; }
2383  else if ( *t2 == -SYMBOL ) { type2 = 2; sgn2 = 1; }
2384  else if ( ABS(t2[*t2-1]) == *t2-ARGHEAD-1 ) {
2385  type2 = 3; if ( t2[*t2-1] < 0 ) sgn2 = -1; else sgn2 = 1; }
2386  else { type2 = 4; if ( t2[*t2-1] < 0 ) sgn2 = -1; else sgn2 = 1; }
2387  if ( type1 > type2 ) {
2388  t3 = t1; t1 = t2; t2 = t3;
2389  type3 = type1; type1 = type2; type2 = type3;
2390  sgn3 = sgn1; sgn1 = sgn2; sgn2 = sgn3;
2391  }
2392  nosymbols = 1; sgn3 = 1;
2393  switch ( type1 ) {
2394  case 1:
2395  if ( type2 == 1 ) {
2396  x = sgn2 * t1[1];
2397  x += sgn1 * t2[1];
2398  if ( x > MAXPOSITIVE || x < -(MAXPOSITIVE+1) ) {
2399  if ( x < 0 ) { sgn1 = -3; x = -x; }
2400  else sgn1 = 3;
2401  *to++ = ARGHEAD+4;
2402  *to++ = 0;
2403  FILLARG(to)
2404  *to++ = 4; *to++ = (UWORD)x; *to++ = 1; *to++ = sgn1;
2405  }
2406  else { *to++ = -SNUMBER; *to++ = (WORD)x; }
2407  }
2408  else if ( type2 == 2 ) {
2409  *to++ = ARGHEAD+8; *to++ = 0; FILLARG(to)
2410  *to++ = 8; *to++ = SYMBOL; *to++ = 4; *to++ = t2[1]; *to++ = 1;
2411  *to++ = ABS(t1[1])+1;
2412  *to++ = 1;
2413  *to++ = 3*sgn1;
2414  }
2415  else if ( type2 == 3 ) {
2416  tt1 = (WORD *)scratch; tt1[0] = ABS(t1[1]); size1 = 1;
2417  tt2 = t2+ARGHEAD+1; size2 = (ABS(t2[*t2-1])-1)/2;
2418  t3 = to;
2419  *to++ = 0; *to++ = 0; FILLARG(to) *to++ = 0;
2420  goto DoCoeffi;
2421  }
2422  else {
2423 /*
2424  t1 is (short) numeric, t2 has the symbol(s).
2425 */
2426  tt1 = (WORD *)scratch; tt1[0] = ABS(t1[1]); size1 = 1;
2427  tt2 = t2+ARGHEAD+1; tt2 += tt2[1]; size2 = (ABS(t2[*t2-1])-1)/2;
2428  t3 = to; i = tt2 - t2; r = t2;
2429  NCOPY(to,r,i)
2430  nosymbols = 0;
2431  goto DoCoeffi;
2432  }
2433  break;
2434  case 2:
2435  if ( type2 == 2 ) {
2436  if ( t1[1] == t2[1] ) {
2437  if ( ( symbols[t1[1]].maxpower == 4 )
2438  && ( ( symbols[t1[1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) ) {
2439  *to++ = -SNUMBER; *to++ = -2;
2440  }
2441  else if ( symbols[t1[1]].maxpower == 2 ) {
2442  *to++ = -SNUMBER; *to++ = 2;
2443  }
2444  else {
2445  *to++ = ARGHEAD+8; *to++ = 0; FILLARG(to)
2446  *to++ = 8; *to++ = SYMBOL; *to++ = 4;
2447  *to++ = t1[1]; *to++ = 2;
2448  *to++ = 2; *to++ = 1; *to++ = 3;
2449  }
2450  }
2451  else {
2452  *to++ = ARGHEAD+10; *to++ = 0; FILLARG(to)
2453  *to++ = 10; *to++ = SYMBOL; *to++ = 6;
2454  if ( t1[1] < t2[1] ) {
2455  *to++ = t1[1]; *to++ = 1; *to++ = t2[1]; *to++ = 1;
2456  }
2457  else {
2458  *to++ = t2[1]; *to++ = 1; *to++ = t1[1]; *to++ = 1;
2459  }
2460  *to++ = 2; *to++ = 1; *to++ = 3;
2461  }
2462  }
2463  else if ( type2 == 3 ) {
2464  t3 = to;
2465  *to++ = 0; *to++ = 0; FILLARG(to) *to++ = 0;
2466  *to++ = SYMBOL; *to++ = 4; *to++ = t1[1]; *to++ = 1;
2467  tt1 = scratch; tt1[1] = 1; size1 = 1;
2468  tt2 = t2+ARGHEAD+1; size2 = (ABS(t2[*t2-1])-1)/2;
2469  nosymbols = 0;
2470  goto DoCoeffi;
2471  }
2472  else {
2473  tt1 = scratch; tt1[0] = 1; size1 = 1;
2474  t3 = to;
2475  *to++ = 0; *to++ = 0; FILLARG(to) *to++ = 0;
2476  *to++ = SYMBOL; *to++ = 0;
2477  tt2 = t2 + ARGHEAD+3; it2 = tt2[-1]-2;
2478  while ( it2 > 0 ) {
2479  if ( *tt2 == t1[1] ) {
2480  pow = tt2[1]+1;
2481  root = symbols[*tt2].maxpower;
2482  if ( pow >= root ) pow -= root;
2483  if ( ( symbols[*tt2].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
2484  if ( ( root & 1 ) == 0 && pow >= root/2 ) {
2485  pow -= root/2; sgn3 = -sgn3;
2486  }
2487  }
2488  if ( pow != 0 ) {
2489  *to++ = *tt2; *to++ = pow;
2490  }
2491  tt2 += 2; it2 -= 2;
2492  break;
2493  }
2494  else if ( t1[1] < *tt2 ) {
2495  *to++ = t1[1]; *to++ = 1; break;
2496  }
2497  else {
2498  *to++ = *tt2++; *to++ = *tt2++; it2 -= 2;
2499  if ( it2 <= 0 ) { *to++ = t1[1]; *to++ = 1; }
2500  }
2501  }
2502  while ( it2 > 0 ) { *to++ = *tt2++; *to++ = *tt2++; it2 -= 2; }
2503  if ( (to - t3) > ARGHEAD+3 ) {
2504  t3[ARGHEAD+2] = (to-t3)-ARGHEAD-1; /* size of the SYMBOL field */
2505  nosymbols = 0;
2506  }
2507  else {
2508  to = t3+ARGHEAD+1; /* no SYMBOL field */
2509  }
2510  size2 = (ABS(t2[*t2-1])-1)/2;
2511  goto DoCoeffi;
2512  }
2513  break;
2514  case 3:
2515  if ( type2 == 3 ) {
2516 /*
2517  Both are numeric
2518 */
2519  tt1 = t1+ARGHEAD+1; size1 = (ABS(t1[*t1-1])-1)/2;
2520  tt2 = t2+ARGHEAD+1; size2 = (ABS(t2[*t2-1])-1)/2;
2521  t3 = to;
2522  *to++ = 0; *to++ = 0; FILLARG(to) *to++ = 0;
2523  goto DoCoeffi;
2524  }
2525  else {
2526 /*
2527  t1 is (long) numeric, t2 has the symbol(s).
2528 */
2529  tt1 = t1+ARGHEAD+1; size1 = (ABS(t1[*t1-1])-1)/2;
2530  tt2 = t2+ARGHEAD+1; tt2 += tt2[1]; size2 = (ABS(t2[*t2-1])-1)/2;
2531  t3 = to; i = tt2 - t2; r = t2;
2532  NCOPY(to,r,i)
2533  nosymbols = 0;
2534  goto DoCoeffi;
2535  }
2536  break;
2537  case 4:
2538 /*
2539  Both have roots of unity
2540  1: Merge the lists and simplify if possible
2541 */
2542  tt1 = t1+ARGHEAD+3; it1 = tt1[-1]-2;
2543  tt2 = t2+ARGHEAD+3; it2 = tt2[-1]-2;
2544  t3 = to;
2545  *to++ = 0; *to++ = 0; FILLARG(to)
2546  *to++ = 0; *to++ = SYMBOL; *to++ = 0;
2547  while ( it1 > 0 && it2 > 0 ) {
2548  if ( *tt1 == *tt2 ) {
2549  pow = tt1[1]+tt2[1];
2550  root = symbols[*tt1].maxpower;
2551  if ( pow >= root ) pow -= root;
2552  if ( ( symbols[*tt1].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
2553  if ( ( root & 1 ) == 0 && pow >= root/2 ) {
2554  pow -= root/2; sgn3 = -sgn3;
2555  }
2556  }
2557  if ( pow != 0 ) {
2558  *to++ = *tt1; *to++ = pow;
2559  }
2560  tt1 += 2; tt2 += 2; it1 -= 2; it2 -= 2;
2561  }
2562  else if ( *tt1 < *tt2 ) {
2563  *to++ = *tt1++; *to++ = *tt1++; it1 -= 2;
2564  }
2565  else {
2566  *to++ = *tt2++; *to++ = *tt2++; it2 -= 2;
2567  }
2568  }
2569  while ( it1 > 0 ) { *to++ = *tt1++; *to++ = *tt1++; it1 -= 2; }
2570  while ( it2 > 0 ) { *to++ = *tt2++; *to++ = *tt2++; it2 -= 2; }
2571  if ( (to - t3) > ARGHEAD+3 ) {
2572  t3[ARGHEAD+2] = (to-t3)-ARGHEAD-1; /* size of the SYMBOL field */
2573  nosymbols = 0;
2574  }
2575  else {
2576  to = t3+ARGHEAD+1; /* no SYMBOL field */
2577  }
2578  size1 = (ABS(t1[*t1-1])-1)/2;
2579  size2 = (ABS(t2[*t2-1])-1)/2;
2580 /*
2581  Now tt1 and tt2 are pointing at their coefficients.
2582  sgn1 is the sign of 1, sgn2 is the sign of 2 and sgn3 is an extra
2583  overall sign.
2584 */
2585 DoCoeffi:
2586  if ( AddLong((UWORD *)tt1,size1,(UWORD *)tt2,size2,(UWORD *)to,&size3) ) {
2587  MLOCK(ErrorMessageLock);
2588  MesPrint("Called from StuffRootAdd");
2589  MUNLOCK(ErrorMessageLock);
2590  Terminate(-1);
2591  }
2592  sgn = sgn1*sgn2*sgn3;
2593  if ( nosymbols && size3 == 1 ) {
2594  if ( (UWORD)(to[0]) <= MAXPOSITIVE && sgn > 0 ) {
2595  sgn1 = to[0];
2596  to = t3; *to++ = -SNUMBER; *to++ = sgn1;
2597  }
2598  else if ( (UWORD)(to[0]) <= (MAXPOSITIVE+1) && sgn < 0 ) {
2599  sgn1 = to[0];
2600  to = t3; *to++ = -SNUMBER; *to++ = -sgn1;
2601  }
2602  else goto genericcoef;
2603  }
2604  else {
2605 genericcoef:
2606  to += size3;
2607  sgn = sgn*(2*size3+1);
2608  *to++ = 1;
2609  while ( size3 > 1 ) { *to++ = 0; size3--; }
2610  *to++ = sgn;
2611  t3[0] = to - t3;
2612  t3[ARGHEAD] = t3[0] - ARGHEAD;
2613  }
2614  break;
2615  }
2616  return(to);
2617 }
2618 
2619 #endif
2620 
2621 /*
2622  #] StuffRootAdd :
2623 */
#define PHEAD
Definition: ftypes.h:56
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:2865
WORD CompCoef(WORD *, WORD *)
Definition: reken.c:3012