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