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