FORM  4.2.1
argument.c
Go to the documentation of this file.
1 
7 /* #[ License : */
8 /*
9  * Copyright (C) 1984-2017 J.A.M. Vermaseren
10  * When using this file you are requested to refer to the publication
11  * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
12  * This is considered a matter of courtesy as the development was paid
13  * for by FOM the Dutch physics granting agency and we would like to
14  * be able to track its scientific use to convince FOM of its value
15  * for the community.
16  *
17  * This file is part of FORM.
18  *
19  * FORM is free software: you can redistribute it and/or modify it under the
20  * terms of the GNU General Public License as published by the Free Software
21  * Foundation, either version 3 of the License, or (at your option) any later
22  * version.
23  *
24  * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
25  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
26  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
27  * details.
28  *
29  * You should have received a copy of the GNU General Public License along
30  * with FORM. If not, see <http://www.gnu.org/licenses/>.
31  */
32 /* #] License : */
33 
34 /*
35  #[ include : argument.c
36 */
37 
38 #include "form3.h"
39 
40 /*
41  #] include :
42  #[ execarg :
43 
44  Executes the subset of statements in an argument environment.
45  The calling routine should be of the type
46  if ( C->lhs[level][0] == TYPEARG ) {
47  if ( execarg(term,level) ) goto GenCall;
48  level = C->lhs[level][2];
49  goto SkipCount;
50  }
51  Note that there will be cases in which extra space is needed.
52  In addition the compare with C->numlhs isn't very fine, because we
53  need to insert a different value (C->lhs[level][2]).
54 */
55 
56 WORD execarg(PHEAD WORD *term, WORD level)
57 {
58  GETBIDENTITY
59  WORD *t, *r, *m, *v;
60  WORD *start, *stop, *rstop, *r1, *r2 = 0, *r3 = 0, *r4, *r5, *r6, *r7, *r8, *r9;
61  WORD *mm, *mstop, *rnext, *rr, *factor, type, ngcd, nq;
62  CBUF *C = cbuf+AM.rbufnum, *CC = cbuf+AT.ebufnum;
63  WORD i, j, k, oldnumlhs = AR.Cnumlhs, count, action = 0, olddefer = AR.DeferFlag;
64  WORD oldnumrhs = CC->numrhs, size, pow, jj;
65  LONG oldcpointer = CC->Pointer - CC->Buffer, oldppointer = AT.pWorkPointer, lp;
66  WORD *oldwork = AT.WorkPointer, *oldwork2, scale, renorm;
67  WORD kLCM = 0, kGCD = 0, kGCD2, kkLCM = 0, jLCM = 0, jGCD, sign = 1;
68  int ii;
69  UWORD *EAscrat, *GCDbuffer = 0, *GCDbuffer2, *LCMbuffer, *LCMb, *LCMc;
70  AT.WorkPointer += *term;
71  start = C->lhs[level];
72  AR.Cnumlhs = start[2];
73  stop = start + start[1];
74  type = *start;
75  scale = start[4];
76  renorm = start[5];
77  start += TYPEARGHEADSIZE;
78 /*
79  #[ Dollars :
80 */
81  if ( renorm && start[1] != 0 ) {/* We have to evaluate $ symbols inside () */
82  t = start+1; factor = oldwork2 = v = AT.WorkPointer;
83  i = *t; t++;
84  *v++ = i+3; i--; NCOPY(v,t,i);
85  *v++ = 1; *v++ = 1; *v++ = 3;
86  AT.WorkPointer = v;
87  start = t; AR.Eside = LHSIDEX;
88  NewSort(BHEAD0);
89  if ( Generator(BHEAD factor,AR.Cnumlhs) ) {
91  AT.WorkPointer = oldwork;
92  return(-1);
93  }
94  AT.WorkPointer = v;
95  if ( EndSort(BHEAD factor,0) < 0 ) {}
96  if ( *factor && *(factor+*factor) != 0 ) {
97  MLOCK(ErrorMessageLock);
98  MesPrint("&$ in () does not evaluate into a single term");
99  MUNLOCK(ErrorMessageLock);
100  return(-1);
101  }
102  AR.Eside = RHSIDE;
103  if ( *factor > 0 ) {
104  v = factor+*factor;
105  v -= ABS(v[-1]);
106  *factor = v-factor;
107  }
108  AT.WorkPointer = v;
109  }
110  else {
111  if ( *start < 0 ) {
112  factor = start + 1;
113  start += -*start;
114  }
115  else factor = 0;
116  }
117 /*
118  #] Dollars :
119 */
120  t = term;
121  r = t + *t;
122  rstop = r - ABS(r[-1]);
123  t++;
124 /*
125  #[ Argument detection : + argument statement
126 */
127  while ( t < rstop ) {
128  if ( *t >= FUNCTION && functions[*t-FUNCTION].spec == 0 ) {
129 /*
130  We have a function. First count the number of arguments.
131  Tensors are excluded.
132 */
133  count = 0;
134  v = t;
135  m = t + FUNHEAD;
136  r = t + t[1];
137  while ( m < r ) {
138  count++;
139  NEXTARG(m)
140  }
141  if ( count <= 0 ) { t += t[1]; continue; }
142 /*
143  Now we take the arguments one by one and test for a match
144 */
145  for ( i = 1; i <= count; i++ ) {
146  m = start;
147  while ( m < stop ) {
148  r = m + m[1];
149  j = *r++;
150  if ( j > 1 ) {
151  while ( --j > 0 ) {
152  if ( *r == i ) goto RightNum;
153  r++;
154  }
155  m = r;
156  continue;
157  }
158 RightNum:
159  if ( m[1] == 2 ) {
160  m += 2;
161  m += *m;
162  goto HaveTodo;
163  }
164  else {
165  r = m + m[1];
166  m += 2;
167  while ( m < r ) {
168  if ( *m == CSET ) {
169  r1 = SetElements + Sets[m[1]].first;
170  r2 = SetElements + Sets[m[1]].last;
171  while ( r1 < r2 ) {
172  if ( *r1++ == *t ) goto HaveTodo;
173  }
174  }
175  else if ( m[1] == *t ) goto HaveTodo;
176  m += 2;
177  }
178  }
179  m += *m;
180  }
181  continue;
182 HaveTodo:
183 /*
184  If we come here we have to do the argument i (first is 1).
185 */
186  sign = 1;
187  action = 1;
188  v[2] |= DIRTYFLAG;
189  r = t + FUNHEAD;
190  j = i;
191  while ( --j > 0 ) { NEXTARG(r) }
192  if ( ( type == TYPESPLITARG ) || ( type == TYPESPLITFIRSTARG )
193  || ( type == TYPESPLITLASTARG ) ) {
194  if ( *t > FUNCTION && *r > 0 ) {
195  WantAddPointers(2);
196  AT.pWorkSpace[AT.pWorkPointer++] = t;
197  AT.pWorkSpace[AT.pWorkPointer++] = r;
198  }
199  continue;
200  }
201  else if ( type == TYPESPLITARG2 ) {
202  if ( *t > FUNCTION && *r > 0 ) {
203  WantAddPointers(2);
204  AT.pWorkSpace[AT.pWorkPointer++] = t;
205  AT.pWorkSpace[AT.pWorkPointer++] = r;
206  }
207  continue;
208  }
209  else if ( type == TYPEFACTARG || type == TYPEFACTARG2 ) {
210  if ( *t > FUNCTION || *t == DENOMINATOR ) {
211  if ( *r > 0 ) {
212  mm = r + ARGHEAD; mstop = r + *r;
213  if ( mm + *mm < mstop ) {
214  WantAddPointers(2);
215  AT.pWorkSpace[AT.pWorkPointer++] = t;
216  AT.pWorkSpace[AT.pWorkPointer++] = r;
217  continue;
218  }
219  if ( *mm == 1+ABS(mstop[-1]) ) continue;
220  if ( mstop[-3] != 1 || mstop[-2] != 1
221  || mstop[-1] != 3 ) {
222  WantAddPointers(2);
223  AT.pWorkSpace[AT.pWorkPointer++] = t;
224  AT.pWorkSpace[AT.pWorkPointer++] = r;
225  continue;
226  }
227  GETSTOP(mm,mstop); mm++;
228  if ( mm + mm[1] < mstop ) {
229  WantAddPointers(2);
230  AT.pWorkSpace[AT.pWorkPointer++] = t;
231  AT.pWorkSpace[AT.pWorkPointer++] = r;
232  continue;
233  }
234  if ( *mm == SYMBOL && ( mm[1] > 4 ||
235  ( mm[3] != 1 && mm[3] != -1 ) ) ) {
236  WantAddPointers(2);
237  AT.pWorkSpace[AT.pWorkPointer++] = t;
238  AT.pWorkSpace[AT.pWorkPointer++] = r;
239  continue;
240  }
241  else if ( *mm == DOTPRODUCT && ( mm[1] > 5 ||
242  ( mm[4] != 1 && mm[4] != -1 ) ) ) {
243  WantAddPointers(2);
244  AT.pWorkSpace[AT.pWorkPointer++] = t;
245  AT.pWorkSpace[AT.pWorkPointer++] = r;
246  continue;
247  }
248  else if ( ( *mm == DELTA || *mm == VECTOR )
249  && mm[1] > 4 ) {
250  WantAddPointers(2);
251  AT.pWorkSpace[AT.pWorkPointer++] = t;
252  AT.pWorkSpace[AT.pWorkPointer++] = r;
253  continue;
254  }
255  }
256  else if ( factor && *factor == 4 && factor[2] == 1 ) {
257  WantAddPointers(2);
258  AT.pWorkSpace[AT.pWorkPointer++] = t;
259  AT.pWorkSpace[AT.pWorkPointer++] = r;
260  continue;
261  }
262  else if ( factor && *factor == 0
263  && ( *r == -SNUMBER && r[1] != 1 ) ) {
264  WantAddPointers(2);
265  AT.pWorkSpace[AT.pWorkPointer++] = t;
266  AT.pWorkSpace[AT.pWorkPointer++] = r;
267  continue;
268  }
269  else if ( *r == -MINVECTOR ) {
270  WantAddPointers(2);
271  AT.pWorkSpace[AT.pWorkPointer++] = t;
272  AT.pWorkSpace[AT.pWorkPointer++] = r;
273  continue;
274  }
275  }
276  continue;
277  }
278  else if ( type == TYPENORM || type == TYPENORM2 || type == TYPENORM3 || type == TYPENORM4 ) {
279  if ( *r < 0 ) {
280  WORD rone;
281  if ( *r == -MINVECTOR ) { rone = -1; *r = -INDEX; }
282  else if ( *r != -SNUMBER || r[1] == 1 || r[1] == 0 ) continue;
283  else { rone = r[1]; r[1] = 1; }
284 /*
285  Now we must multiply the general coefficient by r[1]
286 */
287  if ( scale && ( factor == 0 || *factor ) ) {
288  action = 1;
289  v[2] |= DIRTYFLAG;
290  if ( rone < 0 ) {
291  if ( type == TYPENORM3 ) k = 1;
292  else k = -1;
293  rone = -rone;
294  }
295  else k = 1;
296  r1 = term + *term;
297  size = r1[-1];
298  size = REDLENG(size);
299  if ( scale > 0 ) {
300  for ( jj = 0; jj < scale; jj++ ) {
301  if ( Mully(BHEAD (UWORD *)rstop,&size,(UWORD *)(&rone),k) )
302  goto execargerr;
303  }
304  }
305  else {
306  for ( jj = 0; jj > scale; jj-- ) {
307  if ( Divvy(BHEAD (UWORD *)rstop,&size,(UWORD *)(&rone),k) )
308  goto execargerr;
309  }
310  }
311  size = INCLENG(size);
312  k = size < 0 ? -size: size;
313  rstop[k-1] = size;
314  *term = (WORD)(rstop - term) + k;
315  }
316  continue;
317  }
318 /*
319  Now we have to find a reference term.
320  If factor is defined and *factor != 0 we have to
321  look for the first term that matches the pattern exactly
322  Otherwise the first term plays this role
323  If its coefficient is not one,
324  we must set up a division of the whole argument by
325  this coefficient, and a multiplication of the term
326  when the type is not equal to TYPENORM2.
327  We first multiply the coefficient of the term.
328  Then we set up the division.
329 
330  First find the magic term
331 */
332  if ( type == TYPENORM4 ) {
333 /*
334  For normalizing everything to integers we have to
335  determine for all elements of this argument the LCM of
336  the denominators and the GCD of the numerators.
337 */
338  GCDbuffer = NumberMalloc("execarg");
339  GCDbuffer2 = NumberMalloc("execarg");
340  LCMbuffer = NumberMalloc("execarg");
341  LCMb = NumberMalloc("execarg"); LCMc = NumberMalloc("execarg");
342  r4 = r + *r;
343  r1 = r + ARGHEAD;
344 /*
345  First take the first term to load up the LCM and the GCD
346 */
347  r2 = r1 + *r1;
348  j = r2[-1];
349  if ( j < 0 ) sign = -1;
350  r3 = r2 - ABS(j);
351  k = REDLENG(j);
352  if ( k < 0 ) k = -k;
353  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
354  for ( kGCD = 0; kGCD < k; kGCD++ ) GCDbuffer[kGCD] = r3[kGCD];
355  k = REDLENG(j);
356  if ( k < 0 ) k = -k;
357  r3 += k;
358  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
359  for ( kLCM = 0; kLCM < k; kLCM++ ) LCMbuffer[kLCM] = r3[kLCM];
360  r1 = r2;
361 /*
362  Now go through the rest of the terms in this argument.
363 */
364  while ( r1 < r4 ) {
365  r2 = r1 + *r1;
366  j = r2[-1];
367  r3 = r2 - ABS(j);
368  k = REDLENG(j);
369  if ( k < 0 ) k = -k;
370  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
371  if ( ( ( GCDbuffer[0] == 1 ) && ( kGCD == 1 ) ) ) {
372 /*
373  GCD is already 1
374 */
375  }
376  else if ( ( ( k != 1 ) || ( r3[0] != 1 ) ) ) {
377  if ( GcdLong(BHEAD GCDbuffer,kGCD,(UWORD *)r3,k,GCDbuffer2,&kGCD2) ) {
378  NumberFree(GCDbuffer,"execarg");
379  NumberFree(GCDbuffer2,"execarg");
380  NumberFree(LCMbuffer,"execarg");
381  NumberFree(LCMb,"execarg"); NumberFree(LCMc,"execarg");
382  goto execargerr;
383  }
384  kGCD = kGCD2;
385  for ( ii = 0; ii < kGCD; ii++ ) GCDbuffer[ii] = GCDbuffer2[ii];
386  }
387  else {
388  kGCD = 1; GCDbuffer[0] = 1;
389  }
390  k = REDLENG(j);
391  if ( k < 0 ) k = -k;
392  r3 += k;
393  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
394  if ( ( ( LCMbuffer[0] == 1 ) && ( kLCM == 1 ) ) ) {
395  for ( kLCM = 0; kLCM < k; kLCM++ )
396  LCMbuffer[kLCM] = r3[kLCM];
397  }
398  else if ( ( k != 1 ) || ( r3[0] != 1 ) ) {
399  if ( GcdLong(BHEAD LCMbuffer,kLCM,(UWORD *)r3,k,LCMb,&kkLCM) ) {
400  NumberFree(GCDbuffer,"execarg"); NumberFree(GCDbuffer2,"execarg");
401  NumberFree(LCMbuffer,"execarg"); NumberFree(LCMb,"execarg"); NumberFree(LCMc,"execarg");
402  goto execargerr;
403  }
404  DivLong((UWORD *)r3,k,LCMb,kkLCM,LCMb,&kkLCM,LCMc,&jLCM);
405  MulLong(LCMbuffer,kLCM,LCMb,kkLCM,LCMc,&jLCM);
406  for ( kLCM = 0; kLCM < jLCM; kLCM++ )
407  LCMbuffer[kLCM] = LCMc[kLCM];
408  }
409  else {} /* LCM doesn't change */
410  r1 = r2;
411  }
412 /*
413  Now put the factor together: GCD/LCM
414 */
415  r3 = (WORD *)(GCDbuffer);
416  if ( kGCD == kLCM ) {
417  for ( jGCD = 0; jGCD < kGCD; jGCD++ )
418  r3[jGCD+kGCD] = LCMbuffer[jGCD];
419  k = kGCD;
420  }
421  else if ( kGCD > kLCM ) {
422  for ( jGCD = 0; jGCD < kLCM; jGCD++ )
423  r3[jGCD+kGCD] = LCMbuffer[jGCD];
424  for ( jGCD = kLCM; jGCD < kGCD; jGCD++ )
425  r3[jGCD+kGCD] = 0;
426  k = kGCD;
427  }
428  else {
429  for ( jGCD = kGCD; jGCD < kLCM; jGCD++ )
430  r3[jGCD] = 0;
431  for ( jGCD = 0; jGCD < kLCM; jGCD++ )
432  r3[jGCD+kLCM] = LCMbuffer[jGCD];
433  k = kLCM;
434  }
435 /* NumberFree(GCDbuffer,"execarg"); GCDbuffer = 0; */
436  NumberFree(GCDbuffer2,"execarg");
437  NumberFree(LCMbuffer,"execarg");
438  NumberFree(LCMb,"execarg"); NumberFree(LCMc,"execarg");
439  j = 2*k+1;
440 /*
441  Now we have to correct the overal factor
442 
443  We have a little problem here.
444  r3 is in GCDbuffer and we returned that.
445  At the same time we still use it.
446  This works as long as each worker has its own TermMalloc
447 */
448  if ( scale && ( factor == 0 || *factor > 0 ) )
449  goto ScaledVariety;
450 /*
451  The if was added 28-nov-2012 to give MakeInteger also
452  the (0) option.
453 */
454  if ( scale && ( factor == 0 || *factor ) ) {
455  size = term[*term-1];
456  size = REDLENG(size);
457  if ( MulRat(BHEAD (UWORD *)rstop,size,(UWORD *)r3,k,
458  (UWORD *)rstop,&size) ) goto execargerr;
459  size = INCLENG(size);
460  k = size < 0 ? -size: size;
461  rstop[k-1] = size*sign;
462  *term = (WORD)(rstop - term) + k;
463  }
464  }
465  else {
466  if ( factor && *factor >= 1 ) {
467  r4 = r + *r;
468  r1 = r + ARGHEAD;
469  while ( r1 < r4 ) {
470  r2 = r1 + *r1;
471  r3 = r2 - ABS(r2[-1]);
472  j = r3 - r1;
473  r5 = factor;
474  if ( j != *r5 ) { r1 = r2; continue; }
475  r5++; r6 = r1+1;
476  while ( --j > 0 ) {
477  if ( *r5 != *r6 ) break;
478  r5++; r6++;
479  }
480  if ( j > 0 ) { r1 = r2; continue; }
481  break;
482  }
483  if ( r1 >= r4 ) continue;
484  }
485  else {
486  r1 = r + ARGHEAD;
487  r2 = r1 + *r1;
488  r3 = r2 - ABS(r2[-1]);
489  }
490  if ( *r3 == 1 && r3[1] == 1 ) {
491  if ( r2[-1] == 3 ) continue;
492  if ( r2[-1] == -3 && type == TYPENORM3 ) continue;
493  }
494  action = 1;
495  v[2] |= DIRTYFLAG;
496  j = r2[-1];
497  k = REDLENG(j);
498  if ( j < 0 ) j = -j;
499  if ( type == TYPENORM && scale && ( factor == 0 || *factor ) ) {
500 /*
501  Now we correct the overal factor
502 */
503 ScaledVariety:;
504  size = term[*term-1];
505  size = REDLENG(size);
506  if ( scale > 0 ) {
507  for ( jj = 0; jj < scale; jj++ ) {
508  if ( MulRat(BHEAD (UWORD *)rstop,size,(UWORD *)r3,k,
509  (UWORD *)rstop,&size) ) goto execargerr;
510  }
511  }
512  else {
513  for ( jj = 0; jj > scale; jj-- ) {
514  if ( DivRat(BHEAD (UWORD *)rstop,size,(UWORD *)r3,k,
515  (UWORD *)rstop,&size) ) goto execargerr;
516  }
517  }
518  size = INCLENG(size);
519  k = size < 0 ? -size: size;
520  rstop[k-1] = size*sign;
521  *term = (WORD)(rstop - term) + k;
522  }
523  }
524 /*
525  We generate a statement for adapting all terms in the
526  argument sucessively
527 */
528  r4 = AddRHS(AT.ebufnum,1);
529  while ( (r4+j+12) > CC->Top ) r4 = DoubleCbuffer(AT.ebufnum,r4,3);
530  *r4++ = j+1;
531  i = (j-1)/2; /* was (j-1)*2 ????? 17-oct-2017 */
532  for ( k = 0; k < i; k++ ) *r4++ = r3[i+k];
533  for ( k = 0; k < i; k++ ) *r4++ = r3[k];
534  if ( ( type == TYPENORM3 ) || ( type == TYPENORM4 ) ) *r4++ = j*sign;
535  else *r4++ = r3[j-1];
536  *r4++ = 0;
537  CC->rhs[CC->numrhs+1] = r4;
538  CC->Pointer = r4;
539  AT.mulpat[5] = CC->numrhs;
540  AT.mulpat[7] = AT.ebufnum;
541  }
542  else if ( type == TYPEARGTOEXTRASYMBOL ) {
543  WORD n;
544  if ( r[0] < 0 ) {
545  /* The argument is in the fast notation. */
546  WORD tmp[MaX(9,FUNHEAD+5)];
547  switch ( r[0] ) {
548  case -SNUMBER:
549  if ( r[1] == 0 ) {
550  tmp[0] = 0;
551  }
552  else {
553  tmp[0] = 4;
554  tmp[1] = ABS(r[1]);
555  tmp[2] = 1;
556  tmp[3] = r[1] > 0 ? 3 : -3;
557  tmp[4] = 0;
558  }
559  break;
560  case -SYMBOL:
561  tmp[0] = 8;
562  tmp[1] = SYMBOL;
563  tmp[2] = 4;
564  tmp[3] = r[1];
565  tmp[4] = 1;
566  tmp[5] = 1;
567  tmp[6] = 1;
568  tmp[7] = 3;
569  tmp[8] = 0;
570  break;
571  case -INDEX:
572  case -VECTOR:
573  case -MINVECTOR:
574  tmp[0] = 7;
575  tmp[1] = INDEX;
576  tmp[2] = 3;
577  tmp[3] = r[1];
578  tmp[4] = 1;
579  tmp[5] = 1;
580  tmp[6] = r[0] != -MINVECTOR ? 3 : -3;
581  tmp[7] = 0;
582  break;
583  default:
584  if ( r[0] <= -FUNCTION ) {
585  tmp[0] = FUNHEAD+4;
586  tmp[1] = -r[0];
587  tmp[2] = FUNHEAD;
588  ZeroFillRange(tmp,3,1+FUNHEAD);
589  tmp[FUNHEAD+1] = 1;
590  tmp[FUNHEAD+2] = 1;
591  tmp[FUNHEAD+3] = 3;
592  tmp[FUNHEAD+4] = 0;
593  break;
594  }
595  else {
596  MLOCK(ErrorMessageLock);
597  MesPrint("Unknown fast notation found (TYPEARGTOEXTRASYMBOL)");
598  MUNLOCK(ErrorMessageLock);
599  return(-1);
600  }
601  }
602  n = FindSubexpression(tmp);
603  }
604  else {
605  /*
606  * NOTE: writing to r[r[0]] is legal. As long as we work
607  * in a part of the term, at least the coefficient of
608  * the term must follow.
609  */
610  WORD old_rr0 = r[r[0]];
611  r[r[0]] = 0; /* zero-terminated */
612  n = FindSubexpression(r+ARGHEAD);
613  r[r[0]] = old_rr0;
614  }
615  /* Put the new argument in the work space. */
616  if ( AT.WorkPointer+2 > AT.WorkTop ) {
617  MLOCK(ErrorMessageLock);
618  MesWork();
619  MUNLOCK(ErrorMessageLock);
620  return(-1);
621  }
622  r1 = AT.WorkPointer;
623  if ( scale ) { /* means "tonumber" */
624  r1[0] = -SNUMBER;
625  r1[1] = n;
626  }
627  else {
628  r1[0] = -SYMBOL;
629  r1[1] = MAXVARIABLES-n;
630  }
631  /* We need r2, r3, m and k to shift the data. */
632  r2 = r + (r[0] > 0 ? r[0] : r[0] <= -FUNCTION ? 1 : 2);
633  r3 = r;
634  m = r1+ARGHEAD+2;
635  k = 2;
636  goto do_shift;
637  }
638  r3 = r;
639  AR.DeferFlag = 0;
640  if ( *r > 0 ) {
641  NewSort(BHEAD0);
642  action = 1;
643  r2 = r + *r;
644  r += ARGHEAD;
645  while ( r < r2 ) { /* Sum over the terms */
646  m = AT.WorkPointer;
647  j = *r;
648  while ( --j >= 0 ) *m++ = *r++;
649  r1 = AT.WorkPointer;
650  AT.WorkPointer = m;
651 /*
652  What to do with dummy indices?
653 */
654  if ( type == TYPENORM || type == TYPENORM2 || type == TYPENORM3 || type == TYPENORM4 ) {
655  if ( MultDo(BHEAD r1,AT.mulpat) ) goto execargerr;
656  AT.WorkPointer = r1 + *r1;
657  }
658  if ( Generator(BHEAD r1,level) ) goto execargerr;
659  AT.WorkPointer = r1;
660  }
661  }
662  else {
663  r2 = r + (( *r <= -FUNCTION ) ? 1:2);
664  r1 = AT.WorkPointer;
665  ToGeneral(r,r1,0);
666  m = r1 + ARGHEAD;
667  AT.WorkPointer = r1 + *r1;
668  NewSort(BHEAD0);
669  action = 1;
670 /*
671  What to do with dummy indices?
672 */
673  if ( type == TYPENORM || type == TYPENORM2 || type == TYPENORM3 || type == TYPENORM4 ) {
674  if ( MultDo(BHEAD m,AT.mulpat) ) goto execargerr;
675  AT.WorkPointer = m + *m;
676  }
677  if ( (*m != 0 ) && Generator(BHEAD m,level) ) goto execargerr;
678  AT.WorkPointer = r1;
679  }
680  if ( EndSort(BHEAD AT.WorkPointer+ARGHEAD,1) < 0 ) goto execargerr;
681  AR.DeferFlag = olddefer;
682 /*
683  Now shift the sorted entity over the old argument.
684 */
685  m = AT.WorkPointer+ARGHEAD;
686  while ( *m ) m += *m;
687  k = WORDDIF(m,AT.WorkPointer);
688  *AT.WorkPointer = k;
689  AT.WorkPointer[1] = 0;
690  if ( ToFast(AT.WorkPointer,AT.WorkPointer) ) {
691  if ( *AT.WorkPointer <= -FUNCTION ) k = 1;
692  else k = 2;
693  }
694 do_shift:
695  if ( *r3 > 0 ) j = k - *r3;
696  else if ( *r3 <= -FUNCTION ) j = k - 1;
697  else j = k - 2;
698 
699  t[1] += j;
700  action = 1;
701  v[2] |= DIRTYFLAG;
702  if ( j > 0 ) {
703  r = m + j;
704  while ( m > AT.WorkPointer ) *--r = *--m;
705  AT.WorkPointer = r;
706  m = term + *term;
707  r = m + j;
708  while ( m > r2 ) *--r = *--m;
709  }
710  else if ( j < 0 ) {
711  r = r2 + j;
712  r1 = term + *term;
713  while ( r2 < r1 ) *r++ = *r2++;
714  }
715  r = r3;
716  m = AT.WorkPointer;
717  NCOPY(r,m,k);
718  *term += j;
719  rstop += j;
720  CC->numrhs = oldnumrhs;
721  CC->Pointer = CC->Buffer + oldcpointer;
722  }
723  }
724  t += t[1];
725  }
726 /*
727  #] Argument detection :
728  #[ SplitArg : + varieties
729 */
730  if ( ( type == TYPESPLITARG || type == TYPESPLITARG2
731  || type == TYPESPLITFIRSTARG || type == TYPESPLITLASTARG ) &&
732  AT.pWorkPointer > oldppointer ) {
733  t = term+1;
734  r1 = AT.WorkPointer + 1;
735  lp = oldppointer;
736  while ( t < rstop ) {
737  if ( lp < AT.pWorkPointer && t == AT.pWorkSpace[lp] ) {
738  v = t;
739  m = t + FUNHEAD;
740  r = t + t[1];
741  r2 = r1; while ( t < m ) *r1++ = *t++;
742  while ( m < r ) {
743  t = m;
744  NEXTARG(m)
745  if ( lp >= AT.pWorkPointer || t != AT.pWorkSpace[lp+1] ) {
746  if ( *t > 0 ) t[1] = 0;
747  while ( t < m ) *r1++ = *t++;
748  continue;
749  }
750 /*
751  Now we have a nontrivial argument that should be done.
752 */
753  lp += 2;
754  action = 1;
755  v[2] |= DIRTYFLAG;
756  r3 = t + *t;
757  t += ARGHEAD;
758  if ( type == TYPESPLITFIRSTARG ) {
759  r4 = r1; r5 = t; r7 = oldwork;
760  *r1++ = *t + ARGHEAD;
761  for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
762  j = 0;
763  while ( t < r3 ) {
764  i = *t;
765  if ( j == 0 ) {
766  NCOPY(r7,t,i)
767  j++;
768  }
769  else {
770  NCOPY(r1,t,i)
771  }
772  }
773  *r4 = r1 - r4;
774  if ( j ) {
775  if ( ToFast(r4,r4) ) {
776  r1 = r4;
777  if ( *r1 > -FUNCTION ) r1++;
778  r1++;
779  }
780  r7 = oldwork;
781  while ( --j >= 0 ) {
782  r4 = r1; i = *r7;
783  *r1++ = i+ARGHEAD; *r1++ = 0;
784  FILLARG(r1);
785  NCOPY(r1,r7,i)
786  if ( ToFast(r4,r4) ) {
787  r1 = r4;
788  if ( *r1 > -FUNCTION ) r1++;
789  r1++;
790  }
791  }
792  }
793  t = r3;
794  }
795  else if ( type == TYPESPLITLASTARG ) {
796  r4 = r1; r5 = t; r7 = oldwork;
797  *r1++ = *t + ARGHEAD;
798  for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
799  j = 0;
800  while ( t < r3 ) {
801  i = *t;
802  if ( t+i >= r3 ) {
803  NCOPY(r7,t,i)
804  j++;
805  }
806  else {
807  NCOPY(r1,t,i)
808  }
809  }
810  *r4 = r1 - r4;
811  if ( j ) {
812  if ( ToFast(r4,r4) ) {
813  r1 = r4;
814  if ( *r1 > -FUNCTION ) r1++;
815  r1++;
816  }
817  r7 = oldwork;
818  while ( --j >= 0 ) {
819  r4 = r1; i = *r7;
820  *r1++ = i+ARGHEAD; *r1++ = 0;
821  FILLARG(r1);
822  NCOPY(r1,r7,i)
823  if ( ToFast(r4,r4) ) {
824  r1 = r4;
825  if ( *r1 > -FUNCTION ) r1++;
826  r1++;
827  }
828  }
829  }
830  t = r3;
831  }
832  else if ( factor == 0 || ( type == TYPESPLITARG2 && *factor == 0 ) ) {
833  while ( t < r3 ) {
834  r4 = r1;
835  *r1++ = *t + ARGHEAD;
836  for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
837  i = *t;
838  while ( --i >= 0 ) *r1++ = *t++;
839  if ( ToFast(r4,r4) ) {
840  r1 = r4;
841  if ( *r1 > -FUNCTION ) r1++;
842  r1++;
843  }
844  }
845  }
846  else if ( type == TYPESPLITARG2 ) {
847 /*
848  Here we better put the pattern matcher at work?
849  Remember: there are no wildcards.
850 */
851  WORD *oRepFunList = AN.RepFunList;
852  WORD *oWildMask = AT.WildMask, *oWildValue = AN.WildValue;
853  AN.WildValue = AT.locwildvalue; AT.WildMask = AT.locwildvalue+2;
854  AN.NumWild = 0;
855  r4 = r1; r5 = t; r7 = oldwork;
856  *r1++ = *t + ARGHEAD;
857  for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
858  j = 0;
859  while ( t < r3 ) {
860  AN.UseFindOnly = 0; oldwork2 = AT.WorkPointer;
861  AN.RepFunList = r1;
862  AT.WorkPointer = r1+AN.RepFunNum+2;
863  i = *t;
864  if ( FindRest(BHEAD t,factor) &&
865  ( AN.UsedOtherFind || FindOnce(BHEAD t,factor) ) ) {
866  NCOPY(r7,t,i)
867  j++;
868  }
869  else if ( factor[0] == FUNHEAD+1 && factor[1] >= FUNCTION ) {
870  WORD *rr1 = t+1, *rr2 = t+i;
871  rr2 -= ABS(rr2[-1]);
872  while ( rr1 < rr2 ) {
873  if ( *rr1 == factor[1] ) break;
874  rr1 += rr1[1];
875  }
876  if ( rr1 < rr2 ) {
877  NCOPY(r7,t,i)
878  j++;
879  }
880  else {
881  NCOPY(r1,t,i)
882  }
883  }
884  else {
885  NCOPY(r1,t,i)
886  }
887  AT.WorkPointer = oldwork2;
888  }
889  AN.RepFunList = oRepFunList;
890  *r4 = r1 - r4;
891  if ( j ) {
892  if ( ToFast(r4,r4) ) {
893  r1 = r4;
894  if ( *r1 > -FUNCTION ) r1++;
895  r1++;
896  }
897  r7 = oldwork;
898  while ( --j >= 0 ) {
899  r4 = r1; i = *r7;
900  *r1++ = i+ARGHEAD; *r1++ = 0;
901  FILLARG(r1);
902  NCOPY(r1,r7,i)
903  if ( ToFast(r4,r4) ) {
904  r1 = r4;
905  if ( *r1 > -FUNCTION ) r1++;
906  r1++;
907  }
908  }
909  }
910  t = r3;
911  AT.WildMask = oWildMask; AN.WildValue = oWildValue;
912  }
913  else {
914 /*
915  This code deals with splitting off a single term
916 */
917  r4 = r1; r5 = t;
918  *r1++ = *t + ARGHEAD;
919  for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
920  j = 0;
921  while ( t < r3 ) {
922  r6 = t + *t; r6 -= ABS(r6[-1]);
923  if ( (r6 - t) == *factor ) {
924  k = *factor - 1;
925  for ( ; k > 0; k-- ) {
926  if ( t[k] != factor[k] ) break;
927  }
928  if ( k <= 0 ) {
929  j = r3 - t; t += *t; continue;
930  }
931  }
932  else if ( (r6 - t) == 1 && *factor == 0 ) {
933  j = r3 - t; t += *t; continue;
934  }
935  i = *t;
936  NCOPY(r1,t,i)
937  }
938  *r4 = r1 - r4;
939  if ( j ) {
940  if ( ToFast(r4,r4) ) {
941  r1 = r4;
942  if ( *r1 > -FUNCTION ) r1++;
943  r1++;
944  }
945  t = r3 - j;
946  r4 = r1;
947  *r1++ = *t + ARGHEAD;
948  for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
949  i = *t;
950  while ( --i >= 0 ) *r1++ = *t++;
951  if ( ToFast(r4,r4) ) {
952  r1 = r4;
953  if ( *r1 > -FUNCTION ) r1++;
954  r1++;
955  }
956  }
957  t = r3;
958  }
959  }
960  r2[1] = r1 - r2;
961  }
962  else {
963  r = t + t[1];
964  while ( t < r ) *r1++ = *t++;
965  }
966  }
967  r = term + *term;
968  while ( t < r ) *r1++ = *t++;
969  m = AT.WorkPointer;
970  i = m[0] = r1 - m;
971  t = term;
972  while ( --i >= 0 ) *t++ = *m++;
973  if ( AT.WorkPointer < m ) AT.WorkPointer = m;
974  }
975 /*
976  #] SplitArg :
977  #[ FACTARG :
978 */
979  if ( ( type == TYPEFACTARG || type == TYPEFACTARG2 ) &&
980  AT.pWorkPointer > oldppointer ) {
981  t = term+1;
982  r1 = AT.WorkPointer + 1;
983  lp = oldppointer;
984  while ( t < rstop ) {
985  if ( lp < AT.pWorkPointer && AT.pWorkSpace[lp] == t ) {
986  v = t;
987  m = t + FUNHEAD;
988  r = t + t[1];
989  r2 = r1; while ( t < m ) *r1++ = *t++;
990  while ( m < r ) {
991  rr = t = m;
992  NEXTARG(m)
993  if ( lp >= AT.pWorkPointer || AT.pWorkSpace[lp+1] != t ) {
994  if ( *t > 0 ) t[1] = 0;
995  while ( t < m ) *r1++ = *t++;
996  continue;
997  }
998 /*
999  Now we have a nontrivial argument that should be studied.
1000  Try to find common factors.
1001 */
1002  lp += 2;
1003  if ( *t < 0 ) {
1004  if ( factor && ( *factor == 0 && *t == -SNUMBER ) ) {
1005  *r1++ = *t++;
1006  if ( *t == 0 ) *r1++ = *t++;
1007  else { *r1++ = 1; t++; }
1008  continue;
1009  }
1010  else if ( factor && *factor == 4 && factor[2] == 1 ) {
1011  if ( *t == -SNUMBER ) {
1012  if ( factor[3] < 0 || t[1] >= 0 ) {
1013  while ( t < m ) *r1++ = *t++;
1014  }
1015  else {
1016  *r1++ = -SNUMBER; *r1++ = -1;
1017  *r1++ = *t++; *r1++ = -*t++;
1018  }
1019  }
1020  else {
1021  while ( t < m ) *r1++ = *t++;
1022  *r1++ = -SNUMBER; *r1++ = 1;
1023  }
1024  continue;
1025  }
1026  else if ( *t == -MINVECTOR ) {
1027  *r1++ = -VECTOR; t++; *r1++ = *t++;
1028  *r1++ = -SNUMBER; *r1++ = -1;
1029  *r1++ = -SNUMBER; *r1++ = 1;
1030  continue;
1031  }
1032  }
1033 /*
1034  Now we have a nontrivial argument
1035 */
1036  r3 = t + *t;
1037  t += ARGHEAD; r5 = t; /* Store starting point */
1038  /* We have terms from r5 to r3 */
1039  if ( r5+*r5 == r3 && factor ) { /* One term only */
1040  if ( *factor == 0 ) {
1041  GETSTOP(t,r6);
1042  r9 = r1; *r1++ = 0; *r1++ = 1;
1043  FILLARG(r1);
1044  *r1++ = (r6-t)+3; t++;
1045  while ( t < r6 ) *r1++ = *t++;
1046  *r1++ = 1; *r1++ = 1; *r1++ = 3;
1047  *r9 = r1-r9;
1048  if ( ToFast(r9,r9) ) {
1049  if ( *r9 <= -FUNCTION ) r1 = r9+1;
1050  else r1 = r9+2;
1051  }
1052  t = r3; continue;
1053  }
1054  if ( factor[0] == 4 && factor[2] == 1 ) {
1055  GETSTOP(t,r6);
1056  r7 = r1; *r1++ = (r6-t)+3+ARGHEAD; *r1++ = 0;
1057  FILLARG(r1);
1058  *r1++ = (r6-t)+3; t++;
1059  while ( t < r6 ) *r1++ = *t++;
1060  *r1++ = 1; *r1++ = 1; *r1++ = 3;
1061  if ( ToFast(r7,r7) ) {
1062  if ( *r7 <= -FUNCTION ) r1 = r7+1;
1063  else r1 = r7+2;
1064  }
1065  if ( r3[-1] < 0 && factor[3] > 0 ) {
1066  *r1++ = -SNUMBER; *r1++ = -1;
1067  if ( r3[-1] == -3 && r3[-2] == 1
1068  && ( r3[-3] & MAXPOSITIVE ) == r3[-3] ) {
1069  *r1++ = -SNUMBER; *r1++ = r3[-3];
1070  }
1071  else {
1072  *r1++ = (r3-r6)+1+ARGHEAD;
1073  *r1++ = 0;
1074  FILLARG(r1);
1075  *r1++ = (r3-r6+1);
1076  while ( t < r3 ) *r1++ = *t++;
1077  r1[-1] = -r1[-1];
1078  }
1079  }
1080  else {
1081  if ( ( r3[-1] == -3 || r3[-1] == 3 )
1082  && r3[-2] == 1
1083  && ( r3[-3] & MAXPOSITIVE ) == r3[-3] ) {
1084  *r1++ = -SNUMBER; *r1++ = r3[-3];
1085  if ( r3[-1] < 0 ) r1[-1] = - r1[-1];
1086  }
1087  else {
1088  *r1++ = (r3-r6)+1+ARGHEAD;
1089  *r1++ = 0;
1090  FILLARG(r1);
1091  *r1++ = (r3-r6+1);
1092  while ( t < r3 ) *r1++ = *t++;
1093  }
1094  }
1095  t = r3; continue;
1096  }
1097  }
1098 /*
1099  Now we take the first term and look for its pieces
1100  inside the other terms.
1101 
1102  It is at this point that a more general factorization
1103  routine could take over (allowing for writing the output
1104  properly of course).
1105 */
1106  if ( AC.OldFactArgFlag == NEWFACTARG ) {
1107  if ( factor == 0 ) {
1108  WORD *oldworkpointer2 = AT.WorkPointer;
1109  AT.WorkPointer = r1 + AM.MaxTer+FUNHEAD;
1110  if ( ArgFactorize(BHEAD t-ARGHEAD,r1) < 0 ) {
1111  MesCall("ExecArg");
1112  return(-1);
1113  }
1114  AT.WorkPointer = oldworkpointer2;
1115  t = r3;
1116  while ( *r1 ) { NEXTARG(r1) }
1117  }
1118  else {
1119  rnext = t + *t;
1120  GETSTOP(t,r6);
1121  t++;
1122  t = r5; pow = 1;
1123  while ( t < r3 ) {
1124  t += *t; if ( t[-1] > 0 ) { pow = 0; break; }
1125  }
1126 /*
1127  We have to add here the code for computing the GCD
1128  and to divide it out.
1129 
1130  #[ Numerical factor :
1131 */
1132  t = r5;
1133  EAscrat = (UWORD *)(TermMalloc("execarg"));
1134  if ( t + *t == r3 ) goto onetermnew;
1135  GETSTOP(t,r6);
1136  ngcd = t[t[0]-1];
1137  i = abs(ngcd)-1;
1138  while ( --i >= 0 ) EAscrat[i] = r6[i];
1139  t += *t;
1140  while ( t < r3 ) {
1141  GETSTOP(t,r6);
1142  i = t[t[0]-1];
1143  if ( AccumGCD(BHEAD EAscrat,&ngcd,(UWORD *)r6,i) ) goto execargerr;
1144  if ( ngcd == 3 && EAscrat[0] == 1 && EAscrat[1] == 1 ) break;
1145  t += *t;
1146  }
1147  if ( ngcd != 3 || EAscrat[0] != 1 || EAscrat[1] != 1 ) {
1148  if ( pow ) ngcd = -ngcd;
1149  t = r5; r9 = r1; *r1++ = t[-ARGHEAD]; *r1++ = 1;
1150  FILLARG(r1); ngcd = REDLENG(ngcd);
1151  while ( t < r3 ) {
1152  GETSTOP(t,r6);
1153  r7 = t; r8 = r1;
1154  while ( r7 < r6) *r1++ = *r7++;
1155  t += *t;
1156  i = REDLENG(t[-1]);
1157  if ( DivRat(BHEAD (UWORD *)r6,i,EAscrat,ngcd,(UWORD *)r1,&nq) ) goto execargerr;
1158  nq = INCLENG(nq);
1159  i = ABS(nq)-1;
1160  r1 += i; *r1++ = nq; *r8 = r1-r8;
1161  }
1162  *r9 = r1-r9;
1163  ngcd = INCLENG(ngcd);
1164  i = ABS(ngcd)-1;
1165  if ( factor && *factor == 0 ) {}
1166  else if ( ( factor && factor[0] == 4 && factor[2] == 1
1167  && factor[3] == -3 ) || pow == 0 ) {
1168  r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1169  FILLARG(r1); *r1++ = i+2;
1170  for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1171  *r1++ = ngcd;
1172  if ( ToFast(r9,r9) ) r1 = r9+2;
1173  }
1174  else if ( factor && factor[0] == 4 && factor[2] == 1
1175  && factor[3] > 0 && pow ) {
1176  if ( ngcd < 0 ) ngcd = -ngcd;
1177  *r1++ = -SNUMBER; *r1++ = -1;
1178  r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1179  FILLARG(r1); *r1++ = i+2;
1180  for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1181  *r1++ = ngcd;
1182  if ( ToFast(r9,r9) ) r1 = r9+2;
1183  }
1184  else {
1185  if ( ngcd < 0 ) ngcd = -ngcd;
1186  if ( pow ) { *r1++ = -SNUMBER; *r1++ = -1; }
1187  if ( ngcd != 3 || EAscrat[0] != 1 || EAscrat[1] != 1 ) {
1188  r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1189  FILLARG(r1); *r1++ = i+2;
1190  for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1191  *r1++ = ngcd;
1192  if ( ToFast(r9,r9) ) r1 = r9+2;
1193  }
1194  }
1195  }
1196 /*
1197  #] Numerical factor :
1198 */
1199  else {
1200 onetermnew:;
1201  if ( factor == 0 || *factor > 2 ) {
1202  if ( pow > 0 ) {
1203  *r1++ = -SNUMBER; *r1++ = -1;
1204  t = r5;
1205  while ( t < r3 ) {
1206  t += *t; t[-1] = -t[-1];
1207  }
1208  }
1209  t = rr; *r1++ = *t++; *r1++ = 1; t++;
1210  COPYARG(r1,t);
1211  while ( t < m ) *r1++ = *t++;
1212  }
1213  }
1214  TermFree(EAscrat,"execarg");
1215  }
1216  }
1217  else { /* AC.OldFactArgFlag is ON */
1218  {
1219  WORD *mnext, ncom;
1220  rnext = t + *t;
1221  GETSTOP(t,r6);
1222  t++;
1223  if ( factor == 0 ) {
1224  while ( t < r6 ) {
1225 /*
1226  #[ SYMBOL :
1227 */
1228  if ( *t == SYMBOL ) {
1229  r7 = t; r8 = t + t[1]; t += 2;
1230  while ( t < r8 ) {
1231  pow = t[1];
1232  mm = rnext;
1233  while ( mm < r3 ) {
1234  mnext = mm + *mm;
1235  GETSTOP(mm,mstop); mm++;
1236  while ( mm < mstop ) {
1237  if ( *mm != SYMBOL ) mm += mm[1];
1238  else break;
1239  }
1240  if ( *mm == SYMBOL ) {
1241  mstop = mm + mm[1]; mm += 2;
1242  while ( *mm != *t && mm < mstop ) mm += 2;
1243  if ( mm >= mstop ) pow = 0;
1244  else if ( pow > 0 && mm[1] > 0 ) {
1245  if ( mm[1] < pow ) pow = mm[1];
1246  }
1247  else if ( pow < 0 && mm[1] < 0 ) {
1248  if ( mm[1] > pow ) pow = mm[1];
1249  }
1250  else pow = 0;
1251  }
1252  else pow = 0;
1253  if ( pow == 0 ) break;
1254  mm = mnext;
1255  }
1256  if ( pow == 0 ) { t += 2; continue; }
1257 /*
1258  We have a factor
1259 */
1260  action = 1; i = pow;
1261  if ( i > 0 ) {
1262  while ( --i >= 0 ) {
1263  *r1++ = -SYMBOL;
1264  *r1++ = *t;
1265  }
1266  }
1267  else {
1268  while ( i++ < 0 ) {
1269  *r1++ = 8 + ARGHEAD;
1270  for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1271  *r1++ = 8; *r1++ = SYMBOL;
1272  *r1++ = 4; *r1++ = *t; *r1++ = -1;
1273  *r1++ = 1; *r1++ = 1; *r1++ = 3;
1274  }
1275  }
1276 /*
1277  Now we have to remove the symbols
1278 */
1279  t[1] -= pow;
1280  mm = rnext;
1281  while ( mm < r3 ) {
1282  mnext = mm + *mm;
1283  GETSTOP(mm,mstop); mm++;
1284  while ( mm < mstop ) {
1285  if ( *mm != SYMBOL ) mm += mm[1];
1286  else break;
1287  }
1288  mstop = mm + mm[1]; mm += 2;
1289  while ( mm < mstop && *mm != *t ) mm += 2;
1290  mm[1] -= pow;
1291  mm = mnext;
1292  }
1293  t += 2;
1294  }
1295  }
1296 /*
1297  #] SYMBOL :
1298  #[ DOTPRODUCT :
1299 */
1300  else if ( *t == DOTPRODUCT ) {
1301  r7 = t; r8 = t + t[1]; t += 2;
1302  while ( t < r8 ) {
1303  pow = t[2];
1304  mm = rnext;
1305  while ( mm < r3 ) {
1306  mnext = mm + *mm;
1307  GETSTOP(mm,mstop); mm++;
1308  while ( mm < mstop ) {
1309  if ( *mm != DOTPRODUCT ) mm += mm[1];
1310  else break;
1311  }
1312  if ( *mm == DOTPRODUCT ) {
1313  mstop = mm + mm[1]; mm += 2;
1314  while ( ( *mm != *t || mm[1] != t[1] )
1315  && mm < mstop ) mm += 3;
1316  if ( mm >= mstop ) pow = 0;
1317  else if ( pow > 0 && mm[2] > 0 ) {
1318  if ( mm[2] < pow ) pow = mm[2];
1319  }
1320  else if ( pow < 0 && mm[2] < 0 ) {
1321  if ( mm[2] > pow ) pow = mm[2];
1322  }
1323  else pow = 0;
1324  }
1325  else pow = 0;
1326  if ( pow == 0 ) break;
1327  mm = mnext;
1328  }
1329  if ( pow == 0 ) { t += 3; continue; }
1330 /*
1331  We have a factor
1332 */
1333  action = 1; i = pow;
1334  if ( i > 0 ) {
1335  while ( --i >= 0 ) {
1336  *r1++ = 9 + ARGHEAD;
1337  for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1338  *r1++ = 9; *r1++ = DOTPRODUCT;
1339  *r1++ = 5; *r1++ = *t; *r1++ = t[1]; *r1++ = 1;
1340  *r1++ = 1; *r1++ = 1; *r1++ = 3;
1341  }
1342  }
1343  else {
1344  while ( i++ < 0 ) {
1345  *r1++ = 9 + ARGHEAD;
1346  for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1347  *r1++ = 9; *r1++ = DOTPRODUCT;
1348  *r1++ = 5; *r1++ = *t; *r1++ = t[1]; *r1++ = -1;
1349  *r1++ = 1; *r1++ = 1; *r1++ = 3;
1350  }
1351  }
1352 /*
1353  Now we have to remove the dotproducts
1354 */
1355  t[2] -= pow;
1356  mm = rnext;
1357  while ( mm < r3 ) {
1358  mnext = mm + *mm;
1359  GETSTOP(mm,mstop); mm++;
1360  while ( mm < mstop ) {
1361  if ( *mm != DOTPRODUCT ) mm += mm[1];
1362  else break;
1363  }
1364  mstop = mm + mm[1]; mm += 2;
1365  while ( mm < mstop && ( *mm != *t
1366  || mm[1] != t[1] ) ) mm += 3;
1367  mm[2] -= pow;
1368  mm = mnext;
1369  }
1370  t += 3;
1371  }
1372  }
1373 /*
1374  #] DOTPRODUCT :
1375  #[ DELTA/VECTOR :
1376 */
1377  else if ( *t == DELTA || *t == VECTOR ) {
1378  r7 = t; r8 = t + t[1]; t += 2;
1379  while ( t < r8 ) {
1380  mm = rnext;
1381  pow = 1;
1382  while ( mm < r3 ) {
1383  mnext = mm + *mm;
1384  GETSTOP(mm,mstop); mm++;
1385  while ( mm < mstop ) {
1386  if ( *mm != *r7 ) mm += mm[1];
1387  else break;
1388  }
1389  if ( *mm == *r7 ) {
1390  mstop = mm + mm[1]; mm += 2;
1391  while ( ( *mm != *t || mm[1] != t[1] )
1392  && mm < mstop ) mm += 2;
1393  if ( mm >= mstop ) pow = 0;
1394  }
1395  else pow = 0;
1396  if ( pow == 0 ) break;
1397  mm = mnext;
1398  }
1399  if ( pow == 0 ) { t += 2; continue; }
1400 /*
1401  We have a factor
1402 */
1403  action = 1;
1404  *r1++ = 8 + ARGHEAD;
1405  for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1406  *r1++ = 8; *r1++ = *r7;
1407  *r1++ = 4; *r1++ = *t; *r1++ = t[1];
1408  *r1++ = 1; *r1++ = 1; *r1++ = 3;
1409 /*
1410  Now we have to remove the delta's/vectors
1411 */
1412  mm = rnext;
1413  while ( mm < r3 ) {
1414  mnext = mm + *mm;
1415  GETSTOP(mm,mstop); mm++;
1416  while ( mm < mstop ) {
1417  if ( *mm != *r7 ) mm += mm[1];
1418  else break;
1419  }
1420  mstop = mm + mm[1]; mm += 2;
1421  while ( mm < mstop && (
1422  *mm != *t || mm[1] != t[1] ) ) mm += 2;
1423  *mm = mm[1] = NOINDEX;
1424  mm = mnext;
1425  }
1426  *t = t[1] = NOINDEX;
1427  t += 2;
1428  }
1429  }
1430 /*
1431  #] DELTA/VECTOR :
1432  #[ INDEX :
1433 */
1434  else if ( *t == INDEX ) {
1435  r7 = t; r8 = t + t[1]; t += 2;
1436  while ( t < r8 ) {
1437  mm = rnext;
1438  pow = 1;
1439  while ( mm < r3 ) {
1440  mnext = mm + *mm;
1441  GETSTOP(mm,mstop); mm++;
1442  while ( mm < mstop ) {
1443  if ( *mm != *r7 ) mm += mm[1];
1444  else break;
1445  }
1446  if ( *mm == *r7 ) {
1447  mstop = mm + mm[1]; mm += 2;
1448  while ( *mm != *t
1449  && mm < mstop ) mm++;
1450  if ( mm >= mstop ) pow = 0;
1451  }
1452  else pow = 0;
1453  if ( pow == 0 ) break;
1454  mm = mnext;
1455  }
1456  if ( pow == 0 ) { t++; continue; }
1457 /*
1458  We have a factor
1459 */
1460  action = 1;
1461 /*
1462  The next looks like an error.
1463  We should have here a VECTOR or INDEX like object
1464 
1465  *r1++ = 7 + ARGHEAD;
1466  for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1467  *r1++ = 7; *r1++ = *r7;
1468  *r1++ = 3; *r1++ = *t;
1469  *r1++ = 1; *r1++ = 1; *r1++ = 3;
1470 
1471  Replace this by: (11-apr-2007)
1472 */
1473  if ( *t < 0 ) { *r1++ = -VECTOR; }
1474  else { *r1++ = -INDEX; }
1475  *r1++ = *t;
1476 /*
1477  Now we have to remove the index
1478 */
1479  *t = NOINDEX;
1480  mm = rnext;
1481  while ( mm < r3 ) {
1482  mnext = mm + *mm;
1483  GETSTOP(mm,mstop); mm++;
1484  while ( mm < mstop ) {
1485  if ( *mm != *r7 ) mm += mm[1];
1486  else break;
1487  }
1488  mstop = mm + mm[1]; mm += 2;
1489  while ( mm < mstop &&
1490  *mm != *t ) mm += 1;
1491  *mm = NOINDEX;
1492  mm = mnext;
1493  }
1494  t += 1;
1495  }
1496  }
1497 /*
1498  #] INDEX :
1499  #[ FUNCTION :
1500 */
1501  else if ( *t >= FUNCTION ) {
1502 /*
1503  In the next code we should actually look inside
1504  the DENOMINATOR or EXPONENT for noncommuting objects
1505 */
1506  if ( *t >= FUNCTION &&
1507  functions[*t-FUNCTION].commute == 0 ) ncom = 0;
1508  else ncom = 1;
1509  if ( ncom ) {
1510  mm = r5 + 1;
1511  while ( mm < t && ( *mm == DUMMYFUN
1512  || *mm == DUMMYTEN ) ) mm += mm[1];
1513  if ( mm < t ) { t += t[1]; continue; }
1514  }
1515  mm = rnext; pow = 1;
1516  while ( mm < r3 ) {
1517  mnext = mm + *mm;
1518  GETSTOP(mm,mstop); mm++;
1519  while ( mm < mstop ) {
1520  if ( *mm == *t && mm[1] == t[1] ) {
1521  for ( i = 2; i < t[1]; i++ ) {
1522  if ( mm[i] != t[i] ) break;
1523  }
1524  if ( i >= t[1] )
1525  { mm += mm[1]; goto nextmterm; }
1526  }
1527  if ( ncom && *mm != DUMMYFUN && *mm != DUMMYTEN )
1528  { pow = 0; break; }
1529  mm += mm[1];
1530  }
1531  if ( mm >= mstop ) pow = 0;
1532  if ( pow == 0 ) break;
1533 nextmterm: mm = mnext;
1534  }
1535  if ( pow == 0 ) { t += t[1]; continue; }
1536 /*
1537  Copy the function
1538 */
1539  action = 1;
1540  *r1++ = t[1] + 4 + ARGHEAD;
1541  for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
1542  *r1++ = t[1] + 4;
1543  for ( i = 0; i < t[1]; i++ ) *r1++ = t[i];
1544  *r1++ = 1; *r1++ = 1; *r1++ = 3;
1545 /*
1546  Now we have to take out the functions
1547 */
1548  mm = rnext;
1549  while ( mm < r3 ) {
1550  mnext = mm + *mm;
1551  GETSTOP(mm,mstop); mm++;
1552  while ( mm < mstop ) {
1553  if ( *mm == *t && mm[1] == t[1] ) {
1554  for ( i = 2; i < t[1]; i++ ) {
1555  if ( mm[i] != t[i] ) break;
1556  }
1557  if ( i >= t[1] ) {
1558  if ( functions[*t-FUNCTION].spec > 0 )
1559  *mm = DUMMYTEN;
1560  else
1561  *mm = DUMMYFUN;
1562  mm += mm[1];
1563  goto nextterm;
1564  }
1565  }
1566  mm += mm[1];
1567  }
1568 nextterm: mm = mnext;
1569  }
1570  if ( functions[*t-FUNCTION].spec > 0 )
1571  *t = DUMMYTEN;
1572  else
1573  *t = DUMMYFUN;
1574  action = 1;
1575  v[2] = DIRTYFLAG;
1576  t += t[1];
1577  }
1578 /*
1579  #] FUNCTION :
1580 */
1581  else {
1582  t += t[1];
1583  }
1584  }
1585  }
1586  t = r5; pow = 1;
1587  while ( t < r3 ) {
1588  t += *t; if ( t[-1] > 0 ) { pow = 0; break; }
1589  }
1590 /*
1591  We have to add here the code for computing the GCD
1592  and to divide it out.
1593 */
1594 /*
1595  #[ Numerical factor :
1596 */
1597  t = r5;
1598  EAscrat = (UWORD *)(TermMalloc("execarg"));
1599  if ( t + *t == r3 ) goto oneterm;
1600  GETSTOP(t,r6);
1601  ngcd = t[t[0]-1];
1602  i = abs(ngcd)-1;
1603  while ( --i >= 0 ) EAscrat[i] = r6[i];
1604  t += *t;
1605  while ( t < r3 ) {
1606  GETSTOP(t,r6);
1607  i = t[t[0]-1];
1608  if ( AccumGCD(BHEAD EAscrat,&ngcd,(UWORD *)r6,i) ) goto execargerr;
1609  if ( ngcd == 3 && EAscrat[0] == 1 && EAscrat[1] == 1 ) break;
1610  t += *t;
1611  }
1612  if ( ngcd != 3 || EAscrat[0] != 1 || EAscrat[1] != 1 ) {
1613  if ( pow ) ngcd = -ngcd;
1614  t = r5; r9 = r1; *r1++ = t[-ARGHEAD]; *r1++ = 1;
1615  FILLARG(r1); ngcd = REDLENG(ngcd);
1616  while ( t < r3 ) {
1617  GETSTOP(t,r6);
1618  r7 = t; r8 = r1;
1619  while ( r7 < r6) *r1++ = *r7++;
1620  t += *t;
1621  i = REDLENG(t[-1]);
1622  if ( DivRat(BHEAD (UWORD *)r6,i,EAscrat,ngcd,(UWORD *)r1,&nq) ) goto execargerr;
1623  nq = INCLENG(nq);
1624  i = ABS(nq)-1;
1625  r1 += i; *r1++ = nq; *r8 = r1-r8;
1626  }
1627  *r9 = r1-r9;
1628  ngcd = INCLENG(ngcd);
1629  i = ABS(ngcd)-1;
1630  if ( factor && *factor == 0 ) {}
1631  else if ( ( factor && factor[0] == 4 && factor[2] == 1
1632  && factor[3] == -3 ) || pow == 0 ) {
1633  r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1634  FILLARG(r1); *r1++ = i+2;
1635  for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1636  *r1++ = ngcd;
1637  if ( ToFast(r9,r9) ) r1 = r9+2;
1638  }
1639  else if ( factor && factor[0] == 4 && factor[2] == 1
1640  && factor[3] > 0 && pow ) {
1641  if ( ngcd < 0 ) ngcd = -ngcd;
1642  *r1++ = -SNUMBER; *r1++ = -1;
1643  r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1644  FILLARG(r1); *r1++ = i+2;
1645  for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1646  *r1++ = ngcd;
1647  if ( ToFast(r9,r9) ) r1 = r9+2;
1648  }
1649  else {
1650  if ( ngcd < 0 ) ngcd = -ngcd;
1651  if ( pow ) { *r1++ = -SNUMBER; *r1++ = -1; }
1652  if ( ngcd != 3 || EAscrat[0] != 1 || EAscrat[1] != 1 ) {
1653  r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1654  FILLARG(r1); *r1++ = i+2;
1655  for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1656  *r1++ = ngcd;
1657  if ( ToFast(r9,r9) ) r1 = r9+2;
1658  }
1659  }
1660  }
1661 /*
1662  #] Numerical factor :
1663 */
1664  else {
1665 oneterm:;
1666  if ( factor == 0 || *factor > 2 ) {
1667  if ( pow > 0 ) {
1668  *r1++ = -SNUMBER; *r1++ = -1;
1669  t = r5;
1670  while ( t < r3 ) {
1671  t += *t; t[-1] = -t[-1];
1672  }
1673  }
1674  t = rr; *r1++ = *t++; *r1++ = 1; t++;
1675  COPYARG(r1,t);
1676  while ( t < m ) *r1++ = *t++;
1677  }
1678  }
1679  TermFree(EAscrat,"execarg");
1680  }
1681  } /* AC.OldFactArgFlag */
1682  }
1683  r2[1] = r1 - r2;
1684  action = 1;
1685  v[2] = DIRTYFLAG;
1686  }
1687  else {
1688  r = t + t[1];
1689  while ( t < r ) *r1++ = *t++;
1690  }
1691  }
1692  r = term + *term;
1693  while ( t < r ) *r1++ = *t++;
1694  m = AT.WorkPointer;
1695  i = m[0] = r1 - m;
1696  t = term;
1697  while ( --i >= 0 ) *t++ = *m++;
1698  if ( AT.WorkPointer < t ) AT.WorkPointer = t;
1699  }
1700 /*
1701  #] FACTARG :
1702 */
1703  AR.Cnumlhs = oldnumlhs;
1704  if ( action && Normalize(BHEAD term) ) goto execargerr;
1705  AT.WorkPointer = oldwork;
1706  if ( AT.WorkPointer < term + *term ) AT.WorkPointer = term + *term;
1707  AT.pWorkPointer = oldppointer;
1708  if ( GCDbuffer ) NumberFree(GCDbuffer,"execarg");
1709  return(action);
1710 execargerr:
1711  AT.WorkPointer = oldwork;
1712  AT.pWorkPointer = oldppointer;
1713  MLOCK(ErrorMessageLock);
1714  MesCall("execarg");
1715  MUNLOCK(ErrorMessageLock);
1716  return(-1);
1717 }
1718 
1719 /*
1720  #] execarg :
1721  #[ execterm :
1722 */
1723 
1724 WORD execterm(PHEAD WORD *term, WORD level)
1725 {
1726  GETBIDENTITY
1727  CBUF *C = cbuf+AM.rbufnum;
1728  WORD oldnumlhs = AR.Cnumlhs;
1729  WORD maxisat = C->lhs[level][2];
1730  WORD *buffer1 = 0;
1731  WORD *oldworkpointer = AT.WorkPointer;
1732  WORD *t1, i;
1733  WORD olddeferflag = AR.DeferFlag, tryterm = 0;
1734  AR.DeferFlag = 0;
1735  do {
1736  AR.Cnumlhs = C->lhs[level][3];
1737  NewSort(BHEAD0);
1738 /*
1739  Normally for function arguments we do not use PolyFun/PolyRatFun.
1740  Hence NewSort sets the corresponding variables to zero.
1741  Here we overwrite that.
1742 */
1743  AN.FunSorts[AR.sLevel]->PolyFlag = ( AR.PolyFun != 0 ) ? AR.PolyFunType: 0;
1744  if ( AR.PolyFun == 0 ) { AN.FunSorts[AR.sLevel]->PolyFlag = 0; }
1745  else if ( AR.PolyFunType == 1 ) { AN.FunSorts[AR.sLevel]->PolyFlag = 1; }
1746  else if ( AR.PolyFunType == 2 ) {
1747  if ( AR.PolyFunExp == 2 ) AN.FunSorts[AR.sLevel]->PolyFlag = 1;
1748  else AN.FunSorts[AR.sLevel]->PolyFlag = 2;
1749  }
1750  if ( buffer1 ) {
1751  term = buffer1;
1752  while ( *term ) {
1753  t1 = oldworkpointer;
1754  i = *term; while ( --i >= 0 ) *t1++ = *term++;
1755  AT.WorkPointer = t1;
1756  if ( Generator(BHEAD oldworkpointer,level) ) goto exectermerr;
1757  }
1758  }
1759  else {
1760  if ( Generator(BHEAD term,level) ) goto exectermerr;
1761  }
1762  if ( buffer1 ) {
1763  if ( tryterm ) { TermFree(buffer1,"buffer in sort statement"); tryterm = 0; }
1764  else { M_free((void *)buffer1,"buffer in sort statement"); }
1765  buffer1 = 0;
1766  }
1767  AN.tryterm = 1;
1768  if ( EndSort(BHEAD (WORD *)((VOID *)(&buffer1)),2) < 0 ) goto exectermerr;
1769  tryterm = AN.tryterm; AN.tryterm = 0;
1770  level = AR.Cnumlhs;
1771  } while ( AR.Cnumlhs < maxisat );
1772  AR.Cnumlhs = oldnumlhs;
1773  AR.DeferFlag = olddeferflag;
1774  term = buffer1;
1775  while ( *term ) {
1776  t1 = oldworkpointer;
1777  i = *term; while ( --i >= 0 ) *t1++ = *term++;
1778  AT.WorkPointer = t1;
1779  if ( Generator(BHEAD oldworkpointer,level) ) goto exectermerr;
1780  }
1781  if ( tryterm ) { TermFree(buffer1,"buffer in term statement"); tryterm = 0; }
1782  else { M_free(buffer1,"buffer in term statement"); }
1783  buffer1 = 0;
1784  AT.WorkPointer = oldworkpointer;
1785  return(0);
1786 exectermerr:
1787  AT.WorkPointer = oldworkpointer;
1788  AR.DeferFlag = olddeferflag;
1789  MLOCK(ErrorMessageLock);
1790  MesCall("execterm");
1791  MUNLOCK(ErrorMessageLock);
1792  return(-1);
1793 }
1794 
1795 /*
1796  #] execterm :
1797  #[ ArgumentImplode :
1798 */
1799 
1800 int ArgumentImplode(PHEAD WORD *term, WORD *thelist)
1801 {
1802  GETBIDENTITY
1803  WORD *liststart, *liststop, *inlist;
1804  WORD *w, *t, *tend, *tstop, *tt, *ttstop, *ttt, ncount, i;
1805  int action = 0;
1806  liststop = thelist + thelist[1];
1807  liststart = thelist + 2;
1808  t = term;
1809  tend = t + *t;
1810  tstop = tend - ABS(tend[-1]);
1811  t++;
1812  while ( t < tstop ) {
1813  if ( *t >= FUNCTION ) {
1814  inlist = liststart;
1815  while ( inlist < liststop && *inlist != *t ) inlist += inlist[1];
1816  if ( inlist < liststop ) {
1817  tt = t; ttstop = t + t[1]; w = AT.WorkPointer;
1818  for ( i = 0; i < FUNHEAD; i++ ) *w++ = *tt++;
1819  while ( tt < ttstop ) {
1820  ncount = 0;
1821  if ( *tt == -SNUMBER && tt[1] == 0 ) {
1822  ncount = 1; ttt = tt; tt += 2;
1823  while ( tt < ttstop && *tt == -SNUMBER && tt[1] == 0 ) {
1824  ncount++; tt += 2;
1825  }
1826  }
1827  if ( ncount > 0 ) {
1828  if ( tt < ttstop && *tt == -SNUMBER && ( tt[1] == 1 || tt[1] == -1 ) ) {
1829  *w++ = -SNUMBER;
1830  *w++ = (ncount+1) * tt[1];
1831  tt += 2;
1832  action = 1;
1833  }
1834  else if ( ( tt[0] == tt[ARGHEAD] + ARGHEAD )
1835  && ( ABS(tt[tt[0]-1]) == 3 )
1836  && ( tt[tt[0]-2] == 1 )
1837  && ( tt[tt[0]-3] == 1 ) ) { /* Single term with coef +/- 1 */
1838  i = *tt; NCOPY(w,tt,i)
1839  w[-3] = ncount+1;
1840  action = 1;
1841  }
1842  else if ( *tt == -SYMBOL ) {
1843  *w++ = ARGHEAD+8;
1844  *w++ = 0;
1845  FILLARG(w)
1846  *w++ = 8;
1847  *w++ = SYMBOL;
1848  *w++ = tt[1];
1849  *w++ = 1;
1850  *w++ = ncount+1; *w++ = 1; *w++ = 3;
1851  tt += 2;
1852  action = 1;
1853  }
1854  else if ( *tt <= -FUNCTION ) {
1855  *w++ = ARGHEAD+FUNHEAD+4;
1856  *w++ = 0;
1857  FILLARG(w)
1858  *w++ = -*tt++;
1859  *w++ = FUNHEAD+4;
1860  FILLFUN(w)
1861  *w++ = ncount+1; *w++ = 1; *w++ = 3;
1862  action = 1;
1863  }
1864  else {
1865  while ( ttt < tt ) *w++ = *ttt++;
1866  if ( tt < ttstop && *tt == -SNUMBER ) {
1867  *w++ = *tt++; *w++ = *tt++;
1868  }
1869  }
1870  }
1871  else if ( *tt <= -FUNCTION ) {
1872  *w++ = *tt++;
1873  }
1874  else if ( *tt < 0 ) {
1875  *w++ = *tt++;
1876  *w++ = *tt++;
1877  }
1878  else {
1879  i = *tt; NCOPY(w,tt,i)
1880  }
1881  }
1882  AT.WorkPointer[1] = w - AT.WorkPointer;
1883  while ( tt < tend ) *w++ = *tt++;
1884  ttt = AT.WorkPointer; tt = t;
1885  while ( ttt < w ) *tt++ = *ttt++;
1886  term[0] = tt - term;
1887  AT.WorkPointer = tt;
1888  tend = tt; tstop = tt - ABS(tt[-1]);
1889  }
1890  }
1891  t += t[1];
1892  }
1893  if ( action ) {
1894  if ( Normalize(BHEAD term) ) return(-1);
1895  }
1896  return(0);
1897 }
1898 
1899 /*
1900  #] ArgumentImplode :
1901  #[ ArgumentExplode :
1902 */
1903 
1904 int ArgumentExplode(PHEAD WORD *term, WORD *thelist)
1905 {
1906  GETBIDENTITY
1907  WORD *liststart, *liststop, *inlist, *old;
1908  WORD *w, *t, *tend, *tstop, *tt, *ttstop, *ttt, ncount, i;
1909  int action = 0;
1910  LONG x;
1911  liststop = thelist + thelist[1];
1912  liststart = thelist + 2;
1913  t = term;
1914  tend = t + *t;
1915  tstop = tend - ABS(tend[-1]);
1916  t++;
1917  while ( t < tstop ) {
1918  if ( *t >= FUNCTION ) {
1919  inlist = liststart;
1920  while ( inlist < liststop && *inlist != *t ) inlist += inlist[1];
1921  if ( inlist < liststop ) {
1922  tt = t; ttstop = t + t[1]; w = AT.WorkPointer;
1923  for ( i = 0; i < FUNHEAD; i++ ) *w++ = *tt++;
1924  while ( tt < ttstop ) {
1925  if ( *tt == -SNUMBER && tt[1] != 0 ) {
1926  if ( tt[1] < AM.MaxTer/((WORD)sizeof(WORD)*4)
1927  && tt[1] > -(AM.MaxTer/((WORD)sizeof(WORD)*4))
1928  && ( tt[1] > 1 || tt[1] < -1 ) ) {
1929  ncount = ABS(tt[1]);
1930  while ( ncount > 1 ) {
1931  *w++ = -SNUMBER; *w++ = 0; ncount--;
1932  }
1933  *w++ = -SNUMBER;
1934  if ( tt[1] < 0 ) *w++ = -1;
1935  else *w++ = 1;
1936  tt += 2;
1937  action = 1;
1938  }
1939  else {
1940  *w++ = *tt++; *w++ = *tt++;
1941  }
1942  }
1943  else if ( *tt <= -FUNCTION ) {
1944  *w++ = *tt++;
1945  }
1946  else if ( *tt < 0 ) {
1947  *w++ = *tt++;
1948  *w++ = *tt++;
1949  }
1950  else if ( tt[0] == tt[ARGHEAD]+ARGHEAD ) {
1951  ttt = tt + tt[0] - 1;
1952  i = (ABS(ttt[0])-1)/2;
1953  if ( i > 1 ) {
1954 TooMany: old = AN.currentTerm;
1955  AN.currentTerm = term;
1956  MesPrint("Too many arguments in output of ArgExplode");
1957  MesPrint("Term = %t");
1958  AN.currentTerm = old;
1959  return(-1);
1960  }
1961  if ( ttt[-1] != 1 ) goto NoExplode;
1962  x = ttt[-2];
1963  if ( 2*x > (AT.WorkTop-w)-*term ) goto TooMany;
1964  ncount = x - 1;
1965  while ( ncount > 0 ) {
1966  *w++ = -SNUMBER; *w++ = 0; ncount--;
1967  }
1968  ttt[-2] = 1;
1969  i = *tt; NCOPY(w,tt,i)
1970  action = 1;
1971  }
1972  else {
1973 NoExplode:
1974  i = *tt; NCOPY(w,tt,i)
1975  }
1976  }
1977  AT.WorkPointer[1] = w - AT.WorkPointer;
1978  while ( tt < tend ) *w++ = *tt++;
1979  ttt = AT.WorkPointer; tt = t;
1980  while ( ttt < w ) *tt++ = *ttt++;
1981  term[0] = tt - term;
1982  AT.WorkPointer = tt;
1983  tend = tt; tstop = tt - ABS(tt[-1]);
1984  }
1985  }
1986  t += t[1];
1987  }
1988  if ( action ) {
1989  if ( Normalize(BHEAD term) ) return(-1);
1990  }
1991  return(0);
1992 }
1993 
1994 /*
1995  #] ArgumentExplode :
1996  #[ ArgFactorize :
1997 */
2013 #define NEWORDER
2014 
2015 int ArgFactorize(PHEAD WORD *argin, WORD *argout)
2016 {
2017 /*
2018  #[ step 0 : Declarations and initializations
2019 */
2020  WORD *argfree, *argextra, *argcopy, *t, *tstop, *a, *a1, *a2;
2021 #ifdef NEWORDER
2022  WORD *tt;
2023 #endif
2024  WORD startebuf = cbuf[AT.ebufnum].numrhs,oldword;
2025  WORD oldsorttype = AR.SortType, numargs;
2026  int error = 0, action = 0, i, ii, number, sign = 1;
2027 
2028  *argout = 0;
2029 /*
2030  #] step 0 :
2031  #[ step 1 : Take care of ordering
2032 */
2033  AR.SortType = SORTHIGHFIRST;
2034  if ( oldsorttype != AR.SortType ) {
2035  NewSort(BHEAD0);
2036  oldword = argin[*argin]; argin[*argin] = 0;
2037  t = argin+ARGHEAD;
2038  while ( *t ) {
2039  tstop = t + *t;
2040  if ( AN.ncmod != 0 ) {
2041  if ( AN.ncmod != 1 || ( (WORD)AN.cmod[0] < 0 ) ) {
2042  MLOCK(ErrorMessageLock);
2043  MesPrint("Factorization modulus a number, greater than a WORD not implemented.");
2044  MUNLOCK(ErrorMessageLock);
2045  Terminate(-1);
2046  }
2047  if ( Modulus(t) ) {
2048  MLOCK(ErrorMessageLock);
2049  MesCall("ArgFactorize");
2050  MUNLOCK(ErrorMessageLock);
2051  Terminate(-1);
2052  }
2053  if ( !*t) { t = tstop; continue; }
2054  }
2055  StoreTerm(BHEAD t);
2056  t = tstop;
2057  }
2058  EndSort(BHEAD argin+ARGHEAD,0);
2059  argin[*argin] = oldword;
2060  }
2061 /*
2062  #] step 1 :
2063  #[ step 2 : take out the 'content'.
2064 */
2065  argfree = TakeArgContent(BHEAD argin,argout);
2066  {
2067  a1 = argout;
2068  while ( *a1 ) {
2069  if ( a1[0] == -SNUMBER && ( a1[1] == 1 || a1[1] == -1 ) ) {
2070  if ( a1[1] == -1 ) { sign = -sign; a1[1] = 1; }
2071  if ( a1[2] ) {
2072  a = t = a1+2; while ( *t ) NEXTARG(t);
2073  i = t - a1-2;
2074  t = a1; NCOPY(t,a,i);
2075  *t = 0;
2076  continue;
2077  }
2078  else {
2079  a1[0] = 0;
2080  }
2081  break;
2082  }
2083  else if ( a1[0] == FUNHEAD+ARGHEAD+4 && a1[ARGHEAD] == FUNHEAD+4
2084  && a1[*a1-1] == 3 && a1[*a1-2] == 1 && a1[*a1-3] == 1
2085  && a1[ARGHEAD+1] >= FUNCTION ) {
2086  a = t = a1+*a1; while ( *t ) NEXTARG(t);
2087  i = t - a;
2088  *a1 = -a1[ARGHEAD+1]; t = a1+1; NCOPY(t,a,i);
2089  *t = 0;
2090  }
2091  NEXTARG(a1);
2092  }
2093  }
2094  if ( argfree == 0 ) {
2095  argfree = argin;
2096  }
2097  else if ( argfree[0] == ( argfree[ARGHEAD]+ARGHEAD ) ) {
2098  Normalize(BHEAD argfree+ARGHEAD);
2099  argfree[0] = argfree[ARGHEAD]+ARGHEAD;
2100  argfree[1] = 0;
2101  if ( ( argfree[0] == ARGHEAD+4 ) && ( argfree[ARGHEAD+3] == 3 )
2102  && ( argfree[ARGHEAD+1] == 1 ) && ( argfree[ARGHEAD+2] == 1 ) ) {
2103  goto return0;
2104  }
2105  }
2106  else {
2107 /*
2108  The way we took out objects is rather brutish. We have to
2109  normalize
2110 */
2111  NewSort(BHEAD0);
2112  t = argfree+ARGHEAD;
2113  while ( *t ) {
2114  tstop = t + *t;
2115  Normalize(BHEAD t);
2116  StoreTerm(BHEAD t);
2117  t = tstop;
2118  }
2119  EndSort(BHEAD argfree+ARGHEAD,0);
2120  t = argfree+ARGHEAD;
2121  while ( *t ) t += *t;
2122  *argfree = t - argfree;
2123  }
2124 /*
2125  #] step 2 :
2126  #[ step 3 : look whether we have done this one already.
2127 */
2128  if ( ( number = FindArg(BHEAD argfree) ) != 0 ) {
2129  if ( number > 0 ) t = cbuf[AT.fbufnum].rhs[number-1];
2130  else t = cbuf[AC.ffbufnum].rhs[-number-1];
2131 /*
2132  Now position on the result. Remember we have in the cache:
2133  inputarg,0,outputargs,0
2134  t is currently at inputarg. *inputarg is always positive.
2135  in principle this holds also for the arguments in the output
2136  but we take no risks here (in case of future developments).
2137 */
2138  t += *t; t++;
2139  tstop = t;
2140  ii = 0;
2141  while ( *tstop ) {
2142  if ( *tstop == -SNUMBER && tstop[1] == -1 ) {
2143  sign = -sign; ii += 2;
2144  }
2145  NEXTARG(tstop);
2146  }
2147  a = argout; while ( *a ) NEXTARG(a);
2148 #ifndef NEWORDER
2149  if ( sign == -1 ) { *a++ = -SNUMBER; *a++ = -1; *a = 0; sign = 1; }
2150 #endif
2151  i = tstop - t - ii;
2152  ii = a - argout;
2153  a2 = a; a1 = a + i;
2154  *a1 = 0;
2155  while ( ii > 0 ) { *--a1 = *--a2; ii--; }
2156  a = argout;
2157  while ( *t ) {
2158  if ( *t == -SNUMBER && t[1] == -1 ) { t += 2; }
2159  else { COPY1ARG(a,t) }
2160  }
2161  goto return0;
2162  }
2163 /*
2164  #] step 3 :
2165  #[ step 4 : invoke ConvertToPoly
2166 
2167  We make a copy first in case there are no factors
2168 */
2169  argcopy = TermMalloc("argcopy");
2170  for ( i = 0; i <= *argfree; i++ ) argcopy[i] = argfree[i];
2171 
2172  tstop = argfree + *argfree;
2173  {
2174  WORD sumcommu = 0;
2175  t = argfree + ARGHEAD;
2176  while ( t < tstop ) {
2177  sumcommu += DoesCommu(t);
2178  t += *t;
2179  }
2180  if ( sumcommu > 1 ) {
2181  MesPrint("ERROR: Cannot factorize an argument with more than one noncommuting object");
2182  Terminate(-1);
2183  }
2184  }
2185  t = argfree + ARGHEAD;
2186 
2187  while ( t < tstop ) {
2188  if ( ( t[1] != SYMBOL ) && ( *t != (ABS(t[*t-1])+1) ) ) {
2189  action = 1; break;
2190  }
2191  t += *t;
2192  }
2193  if ( action ) {
2194  t = argfree + ARGHEAD;
2195  argextra = AT.WorkPointer;
2196  NewSort(BHEAD0);
2197  while ( t < tstop ) {
2198  if ( LocalConvertToPoly(BHEAD t,argextra,startebuf,0) < 0 ) {
2199  error = -1;
2200 getout:
2201  AR.SortType = oldsorttype;
2202  TermFree(argcopy,"argcopy");
2203  if ( argfree != argin ) TermFree(argfree,"argfree");
2204  MesCall("ArgFactorize");
2205  Terminate(-1);
2206  return(-1);
2207  }
2208  StoreTerm(BHEAD argextra);
2209  t += *t; argextra += *argextra;
2210  }
2211  if ( EndSort(BHEAD argfree+ARGHEAD,0) ) { error = -2; goto getout; }
2212  t = argfree + ARGHEAD;
2213  while ( *t > 0 ) t += *t;
2214  *argfree = t - argfree;
2215  }
2216 /*
2217  #] step 4 :
2218  #[ step 5 : If not in the tables, we have to do this by hard work.
2219 */
2220 
2221  a = argout;
2222  while ( *a ) NEXTARG(a);
2223  if ( poly_factorize_argument(BHEAD argfree,a) < 0 ) {
2224  MesCall("ArgFactorize");
2225  error = -1;
2226  }
2227 /*
2228  #] step 5 :
2229  #[ step 6 : use now ConvertFromPoly
2230 
2231  Be careful: there should be more than one argument now.
2232 */
2233  if ( error == 0 && action ) {
2234  a1 = a; NEXTARG(a1);
2235  if ( *a1 != 0 ) {
2236  CBUF *C = cbuf+AC.cbufnum;
2237  CBUF *CC = cbuf+AT.ebufnum;
2238  WORD *oldworkpointer = AT.WorkPointer;
2239  WORD *argcopy2 = TermMalloc("argcopy2"), *a1, *a2;
2240  a1 = a; a2 = argcopy2;
2241  while ( *a1 ) {
2242  if ( *a1 < 0 ) {
2243  if ( *a1 > -FUNCTION ) *a2++ = *a1++;
2244  *a2++ = *a1++; *a2 = 0;
2245  continue;
2246  }
2247  t = a1 + ARGHEAD;
2248  tstop = a1 + *a1;
2249  argextra = AT.WorkPointer;
2250  NewSort(BHEAD0);
2251  while ( t < tstop ) {
2252  if ( ConvertFromPoly(BHEAD t,argextra,numxsymbol,CC->numrhs-startebuf+numxsymbol
2253  ,startebuf-numxsymbol,1) <= 0 ) {
2254  TermFree(argcopy2,"argcopy2");
2255  LowerSortLevel();
2256  error = -3;
2257  goto getout;
2258  }
2259  t += *t;
2260  AT.WorkPointer = argextra + *argextra;
2261 /*
2262  ConvertFromPoly leaves terms with subexpressions. Hence:
2263 */
2264  if ( Generator(BHEAD argextra,C->numlhs) ) {
2265  TermFree(argcopy2,"argcopy2");
2266  LowerSortLevel();
2267  error = -4;
2268  goto getout;
2269  }
2270  }
2271  AT.WorkPointer = oldworkpointer;
2272  if ( EndSort(BHEAD a2+ARGHEAD,0) ) { error = -5; goto getout; }
2273  t = a2+ARGHEAD; while ( *t ) t += *t;
2274  *a2 = t - a2; a2[1] = 0; ZEROARG(a2);
2275  ToFast(a2,a2); NEXTARG(a2);
2276  a1 = tstop;
2277  }
2278  i = a2 - argcopy2;
2279  a2 = argcopy2; a1 = a;
2280  NCOPY(a1,a2,i);
2281  *a1 = 0;
2282  TermFree(argcopy2,"argcopy2");
2283 /*
2284  Erase the entries we made temporarily in cbuf[AT.ebufnum]
2285 */
2286  CC->numrhs = startebuf;
2287  }
2288  else { /* no factorization. recover the argument from before step 3. */
2289  for ( i = 0; i <= *argcopy; i++ ) a[i] = argcopy[i];
2290  }
2291  }
2292 /*
2293  #] step 6 :
2294  #[ step 7 : Add this one to the tables.
2295 
2296  Possibly drop some elements in the tables
2297  when they become too full.
2298 */
2299  if ( error == 0 && AN.ncmod == 0 ) {
2300  if ( InsertArg(BHEAD argcopy,a,0) < 0 ) { error = -1; }
2301  }
2302 /*
2303  #] step 7 :
2304  #[ step 8 : Clean up and return.
2305 
2306  Change the order of the arguments in argout and a.
2307  Use argcopy as spare space.
2308 */
2309  ii = a - argout;
2310  for ( i = 0; i < ii; i++ ) argcopy[i] = argout[i];
2311  a1 = a;
2312  while ( *a1 ) {
2313  if ( *a1 == -SNUMBER && a1[1] < 0 ) {
2314  sign = -sign; a1[1] = -a1[1];
2315  if ( a1[1] == 1 ) {
2316  a2 = a1+2; while ( *a2 ) NEXTARG(a2);
2317  i = a2-a1-2; a2 = a1+2;
2318  NCOPY(a1,a2,i);
2319  *a1 = 0;
2320  }
2321  while ( *a1 ) NEXTARG(a1);
2322  break;
2323  }
2324  else {
2325  if ( *a1 > 0 && *a1 == a1[ARGHEAD]+ARGHEAD && a1[*a1-1] < 0 ) {
2326  a1[*a1-1] = -a1[*a1-1]; sign = -sign;
2327  }
2328  if ( *a1 == ARGHEAD+4 && a1[ARGHEAD+1] == 1 && a1[ARGHEAD+2] == 1 ) {
2329  a2 = a1+ARGHEAD+4; while ( *a2 ) NEXTARG(a2);
2330  i = a2-a1-ARGHEAD-4; a2 = a1+ARGHEAD+4;
2331  NCOPY(a1,a2,i);
2332  *a1 = 0;
2333  break;
2334  }
2335  while ( *a1 ) NEXTARG(a1);
2336  break;
2337  }
2338  NEXTARG(a1);
2339  }
2340  i = a1 - a;
2341  a2 = argout;
2342  NCOPY(a2,a,i);
2343  for ( i = 0; i < ii; i++ ) *a2++ = argcopy[i];
2344 #ifndef NEWORDER
2345  if ( sign == -1 ) { *a2++ = -SNUMBER; *a2++ = -1; sign = 1; }
2346 #endif
2347  *a2 = 0;
2348  TermFree(argcopy,"argcopy");
2349 return0:
2350  if ( argfree != argin ) TermFree(argfree,"argfree");
2351  if ( oldsorttype != AR.SortType ) {
2352  AR.SortType = oldsorttype;
2353  a = argout;
2354  while ( *a ) {
2355  if ( *a > 0 ) {
2356  NewSort(BHEAD0);
2357  oldword = a[*a]; a[*a] = 0;
2358  t = a+ARGHEAD;
2359  while ( *t ) {
2360  tstop = t + *t;
2361  StoreTerm(BHEAD t);
2362  t = tstop;
2363  }
2364  EndSort(BHEAD a+ARGHEAD,0);
2365  a[*a] = oldword;
2366  a += *a;
2367  }
2368  else { NEXTARG(a); }
2369  }
2370  }
2371 #ifdef NEWORDER
2372  t = argout; numargs = 0;
2373  while ( *t ) {
2374  tt = t;
2375  NEXTARG(t);
2376  if ( *tt == ABS(t[-1])+1+ARGHEAD && sign == -1 ) { t[-1] = -t[-1]; sign = 1; }
2377  else if ( *tt == -SNUMBER && sign == -1 ) { tt[1] = -tt[1]; sign = 1; }
2378  numargs++;
2379  }
2380  if ( sign == -1 ) {
2381  *t++ = -SNUMBER; *t++ = -1; *t = 0; sign = 1; numargs++;
2382  }
2383 #else
2384 /*
2385  Now we have to sort the arguments
2386  First have the number of 'nontrivial/nonnumerical' arguments
2387  Then make a piece of code like in FullSymmetrize with that number
2388  of arguments to be symmetrized.
2389  Put a function in front
2390  Call the Symmetrize routine
2391 */
2392  t = argout; numargs = 0;
2393  while ( *t && *t != -SNUMBER && ( *t < 0 || ( ABS(t[*t-1]) != *t-1 ) ) ) {
2394  NEXTARG(t);
2395  numargs++;
2396  }
2397 #endif
2398  if ( numargs > 1 ) {
2399  WORD *Lijst;
2400  WORD x[3];
2401  x[0] = argout[-FUNHEAD];
2402  x[1] = argout[-FUNHEAD+1];
2403  x[2] = argout[-FUNHEAD+2];
2404  while ( *t ) { NEXTARG(t); }
2405  argout[-FUNHEAD] = SQRTFUNCTION;
2406  argout[-FUNHEAD+1] = t-argout+FUNHEAD;
2407  argout[-FUNHEAD+2] = 0;
2408  AT.WorkPointer = t+1;
2409  Lijst = AT.WorkPointer;
2410  for ( i = 0; i < numargs; i++ ) Lijst[i] = i;
2411  AT.WorkPointer += numargs;
2412  error = Symmetrize(BHEAD argout-FUNHEAD,Lijst,numargs,1,SYMMETRIC);
2413  AT.WorkPointer = Lijst;
2414  argout[-FUNHEAD] = x[0];
2415  argout[-FUNHEAD+1] = x[1];
2416  argout[-FUNHEAD+2] = x[2];
2417 #ifdef NEWORDER
2418 /*
2419  Now we have to get a potential numerical argument to the first position
2420 */
2421  tstop = argout; while ( *tstop ) { NEXTARG(tstop); }
2422  t = argout; number = 0;
2423  while ( *t ) {
2424  tt = t; NEXTARG(t);
2425  if ( *tt == -SNUMBER ) {
2426  if ( number == 0 ) break;
2427  x[0] = tt[1];
2428  while ( tt > argout ) { *--t = *--tt; }
2429  argout[0] = -SNUMBER; argout[1] = x[0];
2430  break;
2431  }
2432  else if ( *tt == ABS(t[-1])+1+ARGHEAD ) {
2433  if ( number == 0 ) break;
2434  ii = t - tt;
2435  for ( i = 0; i < ii; i++ ) tstop[i] = tt[i];
2436  while ( tt > argout ) { *--t = *--tt; }
2437  for ( i = 0; i < ii; i++ ) argout[i] = tstop[i];
2438  *tstop = 0;
2439  break;
2440  }
2441  number++;
2442  }
2443 #endif
2444  }
2445 /*
2446  #] step 8 :
2447 */
2448  return(error);
2449 }
2450 
2451 /*
2452  #] ArgFactorize :
2453  #[ FindArg :
2454 */
2463 WORD FindArg(PHEAD WORD *a)
2464 {
2465  int number;
2466  if ( AN.ncmod != 0 ) return(0); /* no room for mod stuff */
2467  number = FindTree(AT.fbufnum,a);
2468  if ( number >= 0 ) return(number+1);
2469  number = FindTree(AC.ffbufnum,a);
2470  if ( number >= 0 ) return(-number-1);
2471  return(0);
2472 }
2473 
2474 /*
2475  #] FindArg :
2476  #[ InsertArg :
2477 */
2487 WORD InsertArg(PHEAD WORD *argin, WORD *argout,int par)
2488 {
2489  CBUF *C;
2490  WORD *a, i, bufnum;
2491  if ( par == 0 ) {
2492  bufnum = AT.fbufnum;
2493  C = cbuf+bufnum;
2494  if ( C->numrhs >= (C->maxrhs-2) ) CleanupArgCache(BHEAD AT.fbufnum);
2495  }
2496  else if ( par == 1 ) {
2497  bufnum = AC.ffbufnum;
2498  C = cbuf+bufnum;
2499  }
2500  else { return(-1); }
2501  AddRHS(bufnum,1);
2502  AddNtoC(bufnum,*argin,argin,1);
2503  AddToCB(C,0)
2504  a = argout; while ( *a ) NEXTARG(a);
2505  i = a - argout;
2506  AddNtoC(bufnum,i,argout,2);
2507  AddToCB(C,0)
2508  return(InsTree(bufnum,C->numrhs));
2509 }
2510 
2511 /*
2512  #] InsertArg :
2513  #[ CleanupArgCache :
2514 */
2522 int CleanupArgCache(PHEAD WORD bufnum)
2523 {
2524  CBUF *C = cbuf + bufnum;
2525  COMPTREE *boomlijst = C->boomlijst;
2526  LONG *weights = (LONG *)Malloc1(2*(C->numrhs+1)*sizeof(LONG),"CleanupArgCache");
2527  LONG w, whalf, *extraweights;
2528  WORD *a, *to, *from;
2529  int i,j,k;
2530  for ( i = 1; i <= C->numrhs; i++ ) {
2531  weights[i] = ((LONG)i) * (boomlijst[i].usage);
2532  }
2533 /*
2534  Now sort the weights and determine the halfway weight
2535 */
2536  extraweights = weights+C->numrhs+1;
2537  SortWeights(weights+1,extraweights,C->numrhs);
2538  whalf = weights[C->numrhs/2+1];
2539 /*
2540  We should drop everybody with a weight < whalf.
2541 */
2542  to = C->Buffer;
2543  k = 1;
2544  for ( i = 1; i <= C->numrhs; i++ ) {
2545  from = C->rhs[i]; w = ((LONG)i) * (boomlijst[i].usage);
2546  if ( w >= whalf ) {
2547  if ( i < C->numrhs-1 ) {
2548  if ( to == from ) {
2549  to = C->rhs[i+1];
2550  }
2551  else {
2552  j = C->rhs[i+1] - from;
2553  C->rhs[k] = to;
2554  NCOPY(to,from,j)
2555  }
2556  }
2557  else if ( to == from ) {
2558  to += *to + 1; while ( *to ) NEXTARG(to); to++;
2559  }
2560  else {
2561  a = from; a += *a+1; while ( *a ) NEXTARG(a); a++;
2562  j = a - from;
2563  C->rhs[k] = to;
2564  NCOPY(to,from,j)
2565  }
2566  weights[k++] = boomlijst[i].usage;
2567  }
2568  }
2569  C->numrhs = --k;
2570  C->Pointer = to;
2571 /*
2572  Next we need to rebuild the tree.
2573  Note that this can probably be done much faster by using the
2574  remains of the old tree !!!!!!!!!!!!!!!!
2575 */
2576  ClearTree(AT.fbufnum);
2577  for ( i = 1; i <= k; i++ ) {
2578  InsTree(AT.fbufnum,i);
2579  boomlijst[i].usage = weights[i];
2580  }
2581 /*
2582  And cleanup
2583 */
2584  M_free(weights,"CleanupArgCache");
2585  return(0);
2586 }
2587 
2588 /*
2589  #] CleanupArgCache :
2590  #[ ArgSymbolMerge :
2591 */
2592 
2593 int ArgSymbolMerge(WORD *t1, WORD *t2)
2594 {
2595  WORD *t1e = t1+t1[1];
2596  WORD *t2e = t2+t2[1];
2597  WORD *t1a = t1+2;
2598  WORD *t2a = t2+2;
2599  WORD *t3;
2600  while ( t1a < t1e && t2a < t2e ) {
2601  if ( *t1a < *t2a ) {
2602  if ( t1a[1] >= 0 ) {
2603  t3 = t1a+2;
2604  while ( t3 < t1e ) { t3[-2] = *t3; t3[-1] = t3[1]; t3 += 2; }
2605  t1e -= 2;
2606  }
2607  else t1a += 2;
2608  }
2609  else if ( *t1a > *t2a ) {
2610  if ( t2a[1] >= 0 ) t2a += 2;
2611  else {
2612  t3 = t1e;
2613  while ( t3 > t1a ) { *t3 = t3[-2]; t3[1] = t3[-1]; t3 -= 2; }
2614  *t1a++ = *t2a++;
2615  *t1a++ = *t2a++;
2616  t1e += 2;
2617  }
2618  }
2619  else {
2620  if ( t2a[1] < t1a[1] ) t1a[1] = t2a[1];
2621  t1a += 2; t2a += 2;
2622  }
2623  }
2624  while ( t2a < t2e ) {
2625  if ( t2a[1] < 0 ) {
2626  *t1a++ = *t2a++;
2627  *t1a++ = *t2a++;
2628  }
2629  else t2a += 2;
2630  }
2631  while ( t1a < t1e ) {
2632  if ( t1a[1] >= 0 ) {
2633  t3 = t1a+2;
2634  while ( t3 < t1e ) { t3[-2] = *t3; t3[-1] = t3[1]; t3 += 2; }
2635  t1e -= 2;
2636  }
2637  else t1a += 2;
2638  }
2639  t1[1] = t1a - t1;
2640  return(0);
2641 }
2642 
2643 /*
2644  #] ArgSymbolMerge :
2645  #[ ArgDotproductMerge :
2646 */
2647 
2648 int ArgDotproductMerge(WORD *t1, WORD *t2)
2649 {
2650  WORD *t1e = t1+t1[1];
2651  WORD *t2e = t2+t2[1];
2652  WORD *t1a = t1+2;
2653  WORD *t2a = t2+2;
2654  WORD *t3;
2655  while ( t1a < t1e && t2a < t2e ) {
2656  if ( *t1a < *t2a || ( *t1a == *t2a && t1a[1] < t2a[1] ) ) {
2657  if ( t1a[2] >= 0 ) {
2658  t3 = t1a+3;
2659  while ( t3 < t1e ) { t3[-3] = *t3; t3[-2] = t3[1]; t3[-1] = t3[2]; t3 += 3; }
2660  t1e -= 3;
2661  }
2662  else t1a += 3;
2663  }
2664  else if ( *t1a > *t2a || ( *t1a == *t2a && t1a[1] > t2a[1] ) ) {
2665  if ( t2a[2] >= 0 ) t2a += 3;
2666  else {
2667  t3 = t1e;
2668  while ( t3 > t1a ) { *t3 = t3[-3]; t3[1] = t3[-2]; t3[2] = t3[-1]; t3 -= 3; }
2669  *t1a++ = *t2a++;
2670  *t1a++ = *t2a++;
2671  *t1a++ = *t2a++;
2672  t1e += 3;
2673  }
2674  }
2675  else {
2676  if ( t2a[2] < t1a[2] ) t1a[2] = t2a[2];
2677  t1a += 3; t2a += 3;
2678  }
2679  }
2680  while ( t2a < t2e ) {
2681  if ( t2a[2] < 0 ) {
2682  *t1a++ = *t2a++;
2683  *t1a++ = *t2a++;
2684  *t1a++ = *t2a++;
2685  }
2686  else t2a += 3;
2687  }
2688  while ( t1a < t1e ) {
2689  if ( t1a[2] >= 0 ) {
2690  t3 = t1a+3;
2691  while ( t3 < t1e ) { t3[-3] = *t3; t3[-2] = t3[1]; t3[-1] = t3[2]; t3 += 3; }
2692  t1e -= 3;
2693  }
2694  else t1a += 2;
2695  }
2696  t1[1] = t1a - t1;
2697  return(0);
2698 }
2699 
2700 /*
2701  #] ArgDotproductMerge :
2702  #[ TakeArgContent :
2703 */
2716 WORD *TakeArgContent(PHEAD WORD *argin, WORD *argout)
2717 {
2718  GETBIDENTITY
2719  WORD *t, *rnext, *r1, *r2, *r3, *r5, *r6, *r7, *r8, *r9;
2720  WORD pow, *mm, *mnext, *mstop, *argin2 = argin, *argin3 = argin, *argfree;
2721  WORD ncom;
2722  int j, i, act;
2723  r5 = t = argin + ARGHEAD;
2724  r3 = argin + *argin;
2725  rnext = t + *t;
2726  GETSTOP(t,r6);
2727  r1 = argout;
2728  t++;
2729 /*
2730  First pass: arrange everything but the symbols and dotproducts.
2731  They need separate treatment because we have to take out negative
2732  powers.
2733 */
2734  while ( t < r6 ) {
2735 /*
2736  #[ DELTA/VECTOR :
2737 */
2738  if ( *t == DELTA || *t == VECTOR ) {
2739  r7 = t; r8 = t + t[1]; t += 2;
2740  while ( t < r8 ) {
2741  mm = rnext;
2742  pow = 1;
2743  while ( mm < r3 ) {
2744  mnext = mm + *mm;
2745  GETSTOP(mm,mstop); mm++;
2746  while ( mm < mstop ) {
2747  if ( *mm != *r7 ) mm += mm[1];
2748  else break;
2749  }
2750  if ( *mm == *r7 ) {
2751  mstop = mm + mm[1]; mm += 2;
2752  while ( ( *mm != *t || mm[1] != t[1] )
2753  && mm < mstop ) mm += 2;
2754  if ( mm >= mstop ) pow = 0;
2755  }
2756  else pow = 0;
2757  if ( pow == 0 ) break;
2758  mm = mnext;
2759  }
2760  if ( pow == 0 ) { t += 2; continue; }
2761 /*
2762  We have a factor
2763 */
2764  *r1++ = 8 + ARGHEAD;
2765  for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
2766  *r1++ = 8; *r1++ = *r7;
2767  *r1++ = 4; *r1++ = *t; *r1++ = t[1];
2768  *r1++ = 1; *r1++ = 1; *r1++ = 3;
2769  argout = r1;
2770 /*
2771  Now we have to remove the delta's/vectors
2772 */
2773  mm = rnext;
2774  while ( mm < r3 ) {
2775  mnext = mm + *mm;
2776  GETSTOP(mm,mstop); mm++;
2777  while ( mm < mstop ) {
2778  if ( *mm != *r7 ) mm += mm[1];
2779  else break;
2780  }
2781  mstop = mm + mm[1]; mm += 2;
2782  while ( mm < mstop && (
2783  *mm != *t || mm[1] != t[1] ) ) mm += 2;
2784  *mm = mm[1] = NOINDEX;
2785  mm = mnext;
2786  }
2787  *t = t[1] = NOINDEX;
2788  t += 2;
2789  }
2790  }
2791 /*
2792  #] DELTA/VECTOR :
2793  #[ INDEX :
2794 */
2795  else if ( *t == INDEX ) {
2796  r7 = t; r8 = t + t[1]; t += 2;
2797  while ( t < r8 ) {
2798  mm = rnext;
2799  pow = 1;
2800  while ( mm < r3 ) {
2801  mnext = mm + *mm;
2802  GETSTOP(mm,mstop); mm++;
2803  while ( mm < mstop ) {
2804  if ( *mm != *r7 ) mm += mm[1];
2805  else break;
2806  }
2807  if ( *mm == *r7 ) {
2808  mstop = mm + mm[1]; mm += 2;
2809  while ( *mm != *t
2810  && mm < mstop ) mm++;
2811  if ( mm >= mstop ) pow = 0;
2812  }
2813  else pow = 0;
2814  if ( pow == 0 ) break;
2815  mm = mnext;
2816  }
2817  if ( pow == 0 ) { t++; continue; }
2818 /*
2819  We have a factor
2820 */
2821  if ( *t < 0 ) { *r1++ = -VECTOR; }
2822  else { *r1++ = -INDEX; }
2823  *r1++ = *t;
2824  argout = r1;
2825 /*
2826  Now we have to remove the index
2827 */
2828  *t = NOINDEX;
2829  mm = rnext;
2830  while ( mm < r3 ) {
2831  mnext = mm + *mm;
2832  GETSTOP(mm,mstop); mm++;
2833  while ( mm < mstop ) {
2834  if ( *mm != *r7 ) mm += mm[1];
2835  else break;
2836  }
2837  mstop = mm + mm[1]; mm += 2;
2838  while ( mm < mstop &&
2839  *mm != *t ) mm += 1;
2840  *mm = NOINDEX;
2841  mm = mnext;
2842  }
2843  t += 1;
2844  }
2845  }
2846 /*
2847  #] INDEX :
2848  #[ FUNCTION :
2849 */
2850  else if ( *t >= FUNCTION ) {
2851 /*
2852  In the next code we should actually look inside
2853  the DENOMINATOR or EXPONENT for noncommuting objects
2854 */
2855  if ( *t >= FUNCTION &&
2856  functions[*t-FUNCTION].commute == 0 ) ncom = 0;
2857  else ncom = 1;
2858  if ( ncom ) {
2859  mm = r5 + 1;
2860  while ( mm < t && ( *mm == DUMMYFUN
2861  || *mm == DUMMYTEN ) ) mm += mm[1];
2862  if ( mm < t ) { t += t[1]; continue; }
2863  }
2864  mm = rnext; pow = 1;
2865  while ( mm < r3 ) {
2866  mnext = mm + *mm;
2867  GETSTOP(mm,mstop); mm++;
2868  while ( mm < mstop ) {
2869  if ( *mm == *t && mm[1] == t[1] ) {
2870  for ( i = 2; i < t[1]; i++ ) {
2871  if ( mm[i] != t[i] ) break;
2872  }
2873  if ( i >= t[1] )
2874  { mm += mm[1]; goto nextmterm; }
2875  }
2876  if ( ncom && *mm != DUMMYFUN && *mm != DUMMYTEN )
2877  { pow = 0; break; }
2878  mm += mm[1];
2879  }
2880  if ( mm >= mstop ) pow = 0;
2881  if ( pow == 0 ) break;
2882 nextmterm: mm = mnext;
2883  }
2884  if ( pow == 0 ) { t += t[1]; continue; }
2885 /*
2886  Copy the function
2887 */
2888  *r1++ = t[1] + 4 + ARGHEAD;
2889  for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
2890  *r1++ = t[1] + 4;
2891  for ( i = 0; i < t[1]; i++ ) *r1++ = t[i];
2892  *r1++ = 1; *r1++ = 1; *r1++ = 3;
2893  argout = r1;
2894 /*
2895  Now we have to take out the functions
2896 */
2897  mm = rnext;
2898  while ( mm < r3 ) {
2899  mnext = mm + *mm;
2900  GETSTOP(mm,mstop); mm++;
2901  while ( mm < mstop ) {
2902  if ( *mm == *t && mm[1] == t[1] ) {
2903  for ( i = 2; i < t[1]; i++ ) {
2904  if ( mm[i] != t[i] ) break;
2905  }
2906  if ( i >= t[1] ) {
2907  if ( functions[*t-FUNCTION].spec > 0 )
2908  *mm = DUMMYTEN;
2909  else
2910  *mm = DUMMYFUN;
2911  mm += mm[1];
2912  goto nextterm;
2913  }
2914  }
2915  mm += mm[1];
2916  }
2917 nextterm: mm = mnext;
2918  }
2919  if ( functions[*t-FUNCTION].spec > 0 )
2920  *t = DUMMYTEN;
2921  else
2922  *t = DUMMYFUN;
2923  t += t[1];
2924  }
2925 /*
2926  #] FUNCTION :
2927 */
2928  else {
2929  t += t[1];
2930  }
2931  }
2932 /*
2933  #[ SYMBOL :
2934 
2935  Now collect all symbols. We can use the space after r1 as storage
2936 */
2937  t = argin+ARGHEAD;
2938  rnext = t + *t;
2939  r2 = r1;
2940  while ( t < r3 ) {
2941  GETSTOP(t,r6);
2942  t++;
2943  act = 0;
2944  while ( t < r6 ) {
2945  if ( *t == SYMBOL ) {
2946  act = 1;
2947  i = t[1];
2948  NCOPY(r2,t,i)
2949  }
2950  else { t += t[1]; }
2951  }
2952  if ( act == 0 ) {
2953  *r2++ = SYMBOL; *r2++ = 2;
2954  }
2955  t = rnext; rnext = rnext + *rnext;
2956  }
2957  *r2 = 0;
2958  argin2 = argin;
2959 /*
2960  Now we have a list of all symbols as a sequence of SYMBOL subterms.
2961  Any symbol that is absent in a subterm has power zero.
2962  We now need a list of all minimum powers.
2963  This can be done by subsequent merges.
2964 */
2965  r7 = r1; /* The first object into which we merge. */
2966  r8 = r7 + r7[1]; /* The object that gets merged into r7. */
2967  while ( *r8 ) {
2968  r2 = r8 + r8[1]; /* Next object */
2969  ArgSymbolMerge(r7,r8);
2970  r8 = r2;
2971  }
2972 /*
2973  Now we have to divide by the object in r7 and take it apart as factors.
2974  The division can be simple if there are no negative powers.
2975 */
2976  if ( r7[1] > 2 ) {
2977  r8 = r7+2;
2978  r2 = r7 + r7[1];
2979  act = 0;
2980  pow = 0;
2981  while ( r8 < r2 ) {
2982  if ( r8[1] < 0 ) { act = 1; pow += -r8[1]*(ARGHEAD+8); }
2983  else { pow += 2*r8[1]; }
2984  r8 += 2;
2985  }
2986 /*
2987  The amount of space we need to move r7 is given in pow
2988 */
2989  if ( act == 0 ) { /* this can be done 'in situ' */
2990  t = argin + ARGHEAD;
2991  while ( t < r3 ) {
2992  rnext = t + *t;
2993  GETSTOP(t,r6);
2994  t++;
2995  while ( t < r6 ) {
2996  if ( *t != SYMBOL ) { t += t[1]; continue; }
2997  r8 = r7+2; r9 = t + t[1]; t += 2;
2998  while ( ( t < r9 ) && ( r8 < r2 ) ) {
2999  if ( *t == *r8 ) {
3000  t[1] -= r8[1]; t += 2; r8 += 2;
3001  }
3002  else { /* *t must be < than *r8 !!! */
3003  t += 2;
3004  }
3005  }
3006  t = r9;
3007  }
3008  t = rnext;
3009  }
3010 /*
3011  And now the factors that go to argout.
3012  First we have to move r7 out of the way.
3013 */
3014  r8 = r7+pow; i = r7[1];
3015  while ( --i >= 0 ) r8[i] = r7[i];
3016  r2 += pow;
3017  r8 += 2;
3018  while ( r8 < r2 ) {
3019  for ( i = 0; i < r8[1]; i++ ) { *r1++ = -SYMBOL; *r1++ = *r8; }
3020  r8 += 2;
3021  }
3022  }
3023  else { /* this needs a new location */
3024  argin2 = TermMalloc("TakeArgContent2");
3025 /*
3026  We have to multiply the inverse of r7 into argin
3027  The answer should go to argin2.
3028 */
3029  r5 = argin2; *r5++ = 0; *r5++ = 0; FILLARG(r5);
3030  t = argin+ARGHEAD;
3031  while ( t < r3 ) {
3032  rnext = t + *t;
3033  GETSTOP(t,r6);
3034  r9 = r5;
3035  *r5++ = *t++ + r7[1];
3036  while ( t < r6 ) *r5++ = *t++;
3037  i = r7[1] - 2; r8 = r7+2;
3038  *r5++ = r7[0]; *r5++ = r7[1];
3039  while ( i > 0 ) { *r5++ = *r8++; *r5++ = -*r8++; i -= 2; }
3040  while ( t < rnext ) *r5++ = *t++;
3041  Normalize(BHEAD r9);
3042  r5 = r9 + *r9;
3043  }
3044  *r5 = 0;
3045  *argin2 = r5-argin2;
3046 /*
3047  We may have to sort the terms in argin2.
3048 */
3049  NewSort(BHEAD0);
3050  t = argin2+ARGHEAD;
3051  while ( *t ) {
3052  StoreTerm(BHEAD t);
3053  t += *t;
3054  }
3055  t = argin2+ARGHEAD;
3056  if ( EndSort(BHEAD t,0) < 0 ) goto Irreg;
3057  while ( *t ) t += *t;
3058  *argin2 = t - argin2;
3059  r3 = t;
3060 /*
3061  And now the factors that go to argout.
3062  First we have to move r7 out of the way.
3063 */
3064  r8 = r7+pow; i = r7[1];
3065  while ( --i >= 0 ) r8[i] = r7[i];
3066  r2 += pow;
3067  r8 += 2;
3068  while ( r8 < r2 ) {
3069  if ( r8[1] >= 0 ) {
3070  for ( i = 0; i < r8[1]; i++ ) { *r1++ = -SYMBOL; *r1++ = *r8; }
3071  }
3072  else {
3073  for ( i = 0; i < -r8[1]; i++ ) {
3074  *r1++ = ARGHEAD+8; *r1++ = 0;
3075  FILLARG(r1);
3076  *r1++ = 8; *r1++ = SYMBOL; *r1++ = 4; *r1++ = *r8;
3077  *r1++ = -1; *r1++ = 1; *r1++ = 1; *r1++ = 3;
3078  }
3079  }
3080  r8 += 2;
3081  }
3082  argout = r1;
3083  }
3084  }
3085 /*
3086  #] SYMBOL :
3087  #[ DOTPRODUCT :
3088 
3089  Now collect all dotproducts. We can use the space after r1 as storage
3090 */
3091  t = argin2+ARGHEAD;
3092  rnext = t + *t;
3093  r2 = r1;
3094  while ( t < r3 ) {
3095  GETSTOP(t,r6);
3096  t++;
3097  act = 0;
3098  while ( t < r6 ) {
3099  if ( *t == DOTPRODUCT ) {
3100  act = 1;
3101  i = t[1];
3102  NCOPY(r2,t,i)
3103  }
3104  else { t += t[1]; }
3105  }
3106  if ( act == 0 ) {
3107  *r2++ = DOTPRODUCT; *r2++ = 2;
3108  }
3109  t = rnext; rnext = rnext + *rnext;
3110  }
3111  *r2 = 0;
3112  argin3 = argin2;
3113 /*
3114  Now we have a list of all dotproducts as a sequence of DOTPRODUCT
3115  subterms. Any dotproduct that is absent in a subterm has power zero.
3116  We now need a list of all minimum powers.
3117  This can be done by subsequent merges.
3118 */
3119  r7 = r1; /* The first object into which we merge. */
3120  r8 = r7 + r7[1]; /* The object that gets merged into r7. */
3121  while ( *r8 ) {
3122  r2 = r8 + r8[1]; /* Next object */
3123  ArgDotproductMerge(r7,r8);
3124  r8 = r2;
3125  }
3126 /*
3127  Now we have to divide by the object in r7 and take it apart as factors.
3128  The division can be simple if there are no negative powers.
3129 */
3130  if ( r7[1] > 2 ) {
3131  r8 = r7+2;
3132  r2 = r7 + r7[1];
3133  act = 0;
3134  pow = 0;
3135  while ( r8 < r2 ) {
3136  if ( r8[2] < 0 ) { pow += -r8[2]*(ARGHEAD+9); }
3137  else { pow += r8[2]*(ARGHEAD+9); }
3138  r8 += 3;
3139  }
3140 /*
3141  The amount of space we need to move r7 is given in pow
3142  For dotproducts we always need a new location
3143 */
3144  {
3145  argin3 = TermMalloc("TakeArgContent3");
3146 /*
3147  We have to multiply the inverse of r7 into argin
3148  The answer should go to argin2.
3149 */
3150  r5 = argin3; *r5++ = 0; *r5++ = 0; FILLARG(r5);
3151  t = argin2+ARGHEAD;
3152  while ( t < r3 ) {
3153  rnext = t + *t;
3154  GETSTOP(t,r6);
3155  r9 = r5;
3156  *r5++ = *t++ + r7[1];
3157  while ( t < r6 ) *r5++ = *t++;
3158  i = r7[1] - 2; r8 = r7+2;
3159  *r5++ = r7[0]; *r5++ = r7[1];
3160  while ( i > 0 ) { *r5++ = *r8++; *r5++ = *r8++; *r5++ = -*r8++; i -= 3; }
3161  while ( t < rnext ) *r5++ = *t++;
3162  Normalize(BHEAD r9);
3163  r5 = r9 + *r9;
3164  }
3165  *r5 = 0;
3166  *argin3 = r5-argin3;
3167 /*
3168  We may have to sort the terms in argin3.
3169 */
3170  NewSort(BHEAD0);
3171  t = argin3+ARGHEAD;
3172  while ( *t ) {
3173  StoreTerm(BHEAD t);
3174  t += *t;
3175  }
3176  t = argin3+ARGHEAD;
3177  if ( EndSort(BHEAD t,0) < 0 ) goto Irreg;
3178  while ( *t ) t += *t;
3179  *argin3 = t - argin3;
3180  r3 = t;
3181 /*
3182  And now the factors that go to argout.
3183  First we have to move r7 out of the way.
3184 */
3185  r8 = r7+pow; i = r7[1];
3186  while ( --i >= 0 ) r8[i] = r7[i];
3187  r2 += pow;
3188  r8 += 2;
3189  while ( r8 < r2 ) {
3190  for ( i = ABS(r8[2]); i > 0; i-- ) {
3191  *r1++ = ARGHEAD+9; *r1++ = 0; FILLARG(r1);
3192  *r1++ = 9; *r1++ = DOTPRODUCT; *r1++ = 5; *r1++ = *r8;
3193  *r1++ = r8[1]; *r1++ = r8[2] < 0 ? -1: 1;
3194  *r1++ = 1; *r1++ = 1; *r1++ = 3;
3195  }
3196  r8 += 3;
3197  }
3198  argout = r1;
3199  }
3200  }
3201 /*
3202  #] DOTPRODUCT :
3203 
3204  We have now in argin3 the argument stripped of negative powers and
3205  common factors. The only thing left to deal with is to make the
3206  coefficients integer. For that we have to find the LCM of the denominators
3207  and the GCD of the numerators. And to start with, the sign.
3208  We force the sign of the first term to be positive.
3209 */
3210  t = argin3 + ARGHEAD; pow = 1;
3211  t += *t;
3212  if ( t[-1] < 0 ) {
3213  pow = 0;
3214  t[-1] = -t[-1];
3215  while ( t < r3 ) {
3216  t += *t; t[-1] = -t[-1];
3217  }
3218  }
3219 /*
3220  Now the GCD of the numerators and the LCM of the denominators:
3221 */
3222  argfree = TermMalloc("TakeArgContent1");
3223  if ( AN.cmod != 0 ) {
3224  r1 = MakeMod(BHEAD argin3,r1,argfree);
3225  }
3226  else {
3227  r1 = MakeInteger(BHEAD argin3,r1,argfree);
3228  }
3229  if ( pow == 0 ) {
3230  *r1++ = -SNUMBER;
3231  *r1++ = -1;
3232  }
3233  *r1 = 0;
3234 /*
3235  Cleanup
3236 */
3237  if ( argin3 != argin2 ) TermFree(argin3,"TakeArgContent3");
3238  if ( argin2 != argin ) TermFree(argin2,"TakeArgContent2");
3239  return(argfree);
3240 Irreg:
3241  MesPrint("Irregularity while sorting argument in TakeArgContent");
3242  if ( argin3 != argin2 ) TermFree(argin3,"TakeArgContent3");
3243  if ( argin2 != argin ) TermFree(argin2,"TakeArgContent2");
3244  Terminate(-1);
3245  return(0);
3246 }
3247 
3248 /*
3249  #] TakeArgContent :
3250  #[ MakeInteger :
3251 */
3262 WORD *MakeInteger(PHEAD WORD *argin,WORD *argout,WORD *argfree)
3263 {
3264  UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMb, *LCMc;
3265  WORD *r, *r1, *r2, *r3, *r4, *r5, *rnext, i, k, j;
3266  WORD kGCD, kLCM, kGCD2, kkLCM, jLCM, jGCD;
3267  GCDbuffer = NumberMalloc("MakeInteger");
3268  GCDbuffer2 = NumberMalloc("MakeInteger");
3269  LCMbuffer = NumberMalloc("MakeInteger");
3270  LCMb = NumberMalloc("MakeInteger");
3271  LCMc = NumberMalloc("MakeInteger");
3272  r4 = argin + *argin;
3273  r = argin + ARGHEAD;
3274 /*
3275  First take the first term to load up the LCM and the GCD
3276 */
3277  r2 = r + *r;
3278  j = r2[-1];
3279  r3 = r2 - ABS(j);
3280  k = REDLENG(j);
3281  if ( k < 0 ) k = -k;
3282  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
3283  for ( kGCD = 0; kGCD < k; kGCD++ ) GCDbuffer[kGCD] = r3[kGCD];
3284  k = REDLENG(j);
3285  if ( k < 0 ) k = -k;
3286  r3 += k;
3287  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
3288  for ( kLCM = 0; kLCM < k; kLCM++ ) LCMbuffer[kLCM] = r3[kLCM];
3289  r1 = r2;
3290 /*
3291  Now go through the rest of the terms in this argument.
3292 */
3293  while ( r1 < r4 ) {
3294  r2 = r1 + *r1;
3295  j = r2[-1];
3296  r3 = r2 - ABS(j);
3297  k = REDLENG(j);
3298  if ( k < 0 ) k = -k;
3299  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
3300  if ( ( ( GCDbuffer[0] == 1 ) && ( kGCD == 1 ) ) ) {
3301 /*
3302  GCD is already 1
3303 */
3304  }
3305  else if ( ( ( k != 1 ) || ( r3[0] != 1 ) ) ) {
3306  if ( GcdLong(BHEAD GCDbuffer,kGCD,(UWORD *)r3,k,GCDbuffer2,&kGCD2) ) {
3307  NumberFree(GCDbuffer,"MakeInteger");
3308  NumberFree(GCDbuffer2,"MakeInteger");
3309  NumberFree(LCMbuffer,"MakeInteger");
3310  NumberFree(LCMb,"MakeInteger"); NumberFree(LCMc,"MakeInteger");
3311  goto MakeIntegerErr;
3312  }
3313  kGCD = kGCD2;
3314  for ( i = 0; i < kGCD; i++ ) GCDbuffer[i] = GCDbuffer2[i];
3315  }
3316  else {
3317  kGCD = 1; GCDbuffer[0] = 1;
3318  }
3319  k = REDLENG(j);
3320  if ( k < 0 ) k = -k;
3321  r3 += k;
3322  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
3323  if ( ( ( LCMbuffer[0] == 1 ) && ( kLCM == 1 ) ) ) {
3324  for ( kLCM = 0; kLCM < k; kLCM++ )
3325  LCMbuffer[kLCM] = r3[kLCM];
3326  }
3327  else if ( ( k != 1 ) || ( r3[0] != 1 ) ) {
3328  if ( GcdLong(BHEAD LCMbuffer,kLCM,(UWORD *)r3,k,LCMb,&kkLCM) ) {
3329  NumberFree(GCDbuffer,"MakeInteger"); NumberFree(GCDbuffer2,"MakeInteger");
3330  NumberFree(LCMbuffer,"MakeInteger"); NumberFree(LCMb,"MakeInteger"); NumberFree(LCMc,"MakeInteger");
3331  goto MakeIntegerErr;
3332  }
3333  DivLong((UWORD *)r3,k,LCMb,kkLCM,LCMb,&kkLCM,LCMc,&jLCM);
3334  MulLong(LCMbuffer,kLCM,LCMb,kkLCM,LCMc,&jLCM);
3335  for ( kLCM = 0; kLCM < jLCM; kLCM++ )
3336  LCMbuffer[kLCM] = LCMc[kLCM];
3337  }
3338  else {} /* LCM doesn't change */
3339  r1 = r2;
3340  }
3341 /*
3342  Now put the factor together: GCD/LCM
3343 */
3344  r3 = (WORD *)(GCDbuffer);
3345  if ( kGCD == kLCM ) {
3346  for ( jGCD = 0; jGCD < kGCD; jGCD++ )
3347  r3[jGCD+kGCD] = LCMbuffer[jGCD];
3348  k = kGCD;
3349  }
3350  else if ( kGCD > kLCM ) {
3351  for ( jGCD = 0; jGCD < kLCM; jGCD++ )
3352  r3[jGCD+kGCD] = LCMbuffer[jGCD];
3353  for ( jGCD = kLCM; jGCD < kGCD; jGCD++ )
3354  r3[jGCD+kGCD] = 0;
3355  k = kGCD;
3356  }
3357  else {
3358  for ( jGCD = kGCD; jGCD < kLCM; jGCD++ )
3359  r3[jGCD] = 0;
3360  for ( jGCD = 0; jGCD < kLCM; jGCD++ )
3361  r3[jGCD+kLCM] = LCMbuffer[jGCD];
3362  k = kLCM;
3363  }
3364  j = 2*k+1;
3365 /*
3366  Now we have to write this to argout
3367 */
3368  if ( ( j == 3 ) && ( r3[1] == 1 ) && ( (WORD)(r3[0]) > 0 ) ) {
3369  *argout = -SNUMBER;
3370  argout[1] = r3[0];
3371  r1 = argout+2;
3372  }
3373  else {
3374  r1 = argout;
3375  *r1++ = j+1+ARGHEAD; *r1++ = 0; FILLARG(r1);
3376  *r1++ = j+1; r2 = r3;
3377  for ( i = 0; i < k; i++ ) { *r1++ = *r2++; *r1++ = *r2++; }
3378  *r1++ = j;
3379  }
3380 /*
3381  Next we have to take the factor out from the argument.
3382  This cannot be done in location, because the denominator stuff can make
3383  coefficients longer.
3384 */
3385  r2 = argfree + 2; FILLARG(r2)
3386  while ( r < r4 ) {
3387  rnext = r + *r;
3388  j = ABS(rnext[-1]);
3389  r5 = rnext - j;
3390  r3 = r2;
3391  while ( r < r5 ) *r2++ = *r++;
3392  j = (j-1)/2; /* reduced length. Remember, k is the other red length */
3393  if ( DivRat(BHEAD (UWORD *)r5,j,GCDbuffer,k,(UWORD *)r2,&i) ) {
3394  goto MakeIntegerErr;
3395  }
3396  i = 2*i+1;
3397  r2 = r2 + i;
3398  if ( rnext[-1] < 0 ) r2[-1] = -i;
3399  else r2[-1] = i;
3400  *r3 = r2-r3;
3401  r = rnext;
3402  }
3403  *r2 = 0;
3404  argfree[0] = r2-argfree;
3405  argfree[1] = 0;
3406 /*
3407  Cleanup
3408 */
3409  NumberFree(LCMc,"MakeInteger");
3410  NumberFree(LCMb,"MakeInteger");
3411  NumberFree(LCMbuffer,"MakeInteger");
3412  NumberFree(GCDbuffer2,"MakeInteger");
3413  NumberFree(GCDbuffer,"MakeInteger");
3414  return(r1);
3415 
3416 MakeIntegerErr:
3417  MesCall("MakeInteger");
3418  Terminate(-1);
3419  return(0);
3420 }
3421 
3422 /*
3423  #] MakeInteger :
3424  #[ MakeMod :
3425 */
3433 WORD *MakeMod(PHEAD WORD *argin,WORD *argout,WORD *argfree)
3434 {
3435  WORD *r, *instop, *r1, *m, x, xx, ix, ip;
3436  int i;
3437  r = argin; instop = r + *r; r += ARGHEAD;
3438  x = r[*r-3];
3439  if ( r[*r-1] < 0 ) x += AN.cmod[0];
3440  if ( GetModInverses(x,(WORD)(AN.cmod[0]),&ix,&ip) ) {
3441  Terminate(-1);
3442  }
3443  argout[0] = -SNUMBER;
3444  argout[1] = x;
3445  argout[2] = 0;
3446  r1 = argout+2;
3447 /*
3448  Now we have to multiply all coefficients by ix.
3449  This does not make things longer, but we should keep to the conventions
3450  of MakeInteger.
3451 */
3452  m = argfree + ARGHEAD;
3453  while ( r < instop ) {
3454  xx = r[*r-3]; if ( r[*r-1] < 0 ) xx += AN.cmod[0];
3455  xx = (WORD)((((LONG)xx)*ix) % AN.cmod[0]);
3456  if ( xx != 0 ) {
3457  i = *r; NCOPY(m,r,i);
3458  m[-3] = xx; m[-1] = 3;
3459  }
3460  else { r += *r; }
3461  }
3462  *m = 0;
3463  *argfree = m - argfree;
3464  argfree[1] = 0;
3465  argfree += 2; FILLARG(argfree);
3466  return(r1);
3467 }
3468 
3469 /*
3470  #] MakeMod :
3471  #[ SortWeights :
3472 */
3478 void SortWeights(LONG *weights,LONG *extraspace,WORD number)
3479 {
3480  LONG w, *fill, *from1, *from2;
3481  int n1,n2,i;
3482  if ( number >= 4 ) {
3483  n1 = number/2; n2 = number - n1;
3484  SortWeights(weights,extraspace,n1);
3485  SortWeights(weights+n1,extraspace,n2);
3486 /*
3487  We copy the first patch to the extra space. Then we merge
3488  Note that a potential remaining n2 objects are already in place.
3489 */
3490  for ( i = 0; i < n1; i++ ) extraspace[i] = weights[i];
3491  fill = weights; from1 = extraspace; from2 = weights+n1;
3492  while ( n1 > 0 && n2 > 0 ) {
3493  if ( *from1 <= *from2 ) { *fill++ = *from1++; n1--; }
3494  else { *fill++ = *from2++; n2--; }
3495  }
3496  while ( n1 > 0 ) { *fill++ = *from1++; n1--; }
3497  }
3498 /*
3499  Special cases
3500 */
3501  else if ( number == 3 ) { /* 6 permutations of which one is trivial */
3502  if ( weights[0] > weights[1] ) {
3503  if ( weights[1] > weights[2] ) {
3504  w = weights[0]; weights[0] = weights[2]; weights[2] = w;
3505  }
3506  else if ( weights[0] > weights[2] ) {
3507  w = weights[0]; weights[0] = weights[1];
3508  weights[1] = weights[2]; weights[2] = w;
3509  }
3510  else {
3511  w = weights[0]; weights[0] = weights[1]; weights[1] = w;
3512  }
3513  }
3514  else if ( weights[0] > weights[2] ) {
3515  w = weights[0]; weights[0] = weights[2];
3516  weights[2] = weights[1]; weights[1] = w;
3517  }
3518  else if ( weights[1] > weights[2] ) {
3519  w = weights[1]; weights[1] = weights[2]; weights[2] = w;
3520  }
3521  }
3522  else if ( number == 2 ) {
3523  if ( weights[0] > weights[1] ) {
3524  w = weights[0]; weights[0] = weights[1]; weights[1] = w;
3525  }
3526  }
3527  return;
3528 }
3529 
3530 /*
3531  #] SortWeights :
3532 */
int CleanupArgCache(PHEAD WORD bufnum)
Definition: argument.c:2522
WORD * MakeInteger(PHEAD WORD *argin, WORD *argout, WORD *argfree)
Definition: argument.c:3262
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
Definition: notation.c:510
WORD * DoubleCbuffer(int num, WORD *w, int par)
Definition: comtool.c:143
WORD ** lhs
Definition: structs.h:942
Definition: structs.h:938
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition: reken.c:1466
Definition: structs.h:293
WORD * Pointer
Definition: structs.h:941
WORD StoreTerm(PHEAD WORD *)
Definition: sort.c:4332
int usage
Definition: structs.h:299
WORD * TakeArgContent(PHEAD WORD *argin, WORD *argout)
Definition: argument.c:2716
int AddNtoC(int bufnum, int n, WORD *array, int par)
Definition: comtool.c:317
WORD ** rhs
Definition: structs.h:943
void SortWeights(LONG *weights, LONG *extraspace, WORD number)
Definition: argument.c:3478
COMPTREE * boomlijst
Definition: structs.h:948
VOID LowerSortLevel()
Definition: sort.c:4726
#define ZeroFillRange(w, begin, end)
Definition: declare.h:174
WORD InsertArg(PHEAD WORD *argin, WORD *argout, int par)
Definition: argument.c:2487
WORD FindArg(PHEAD WORD *a)
Definition: argument.c:2463
WORD * Buffer
Definition: structs.h:939
WORD NewSort(PHEAD0)
Definition: sort.c:591
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:3072
WORD * MakeMod(PHEAD WORD *argin, WORD *argout, WORD *argfree)
Definition: argument.c:3433
LONG EndSort(PHEAD WORD *, int)
Definition: sort.c:681
WORD * AddRHS(int num, int type)
Definition: comtool.c:214