FORM  4.1
ratio.c
Go to the documentation of this file.
1 
8 /* #[ License : */
9 /*
10  * Copyright (C) 1984-2013 J.A.M. Vermaseren
11  * When using this file you are requested to refer to the publication
12  * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
13  * This is considered a matter of courtesy as the development was paid
14  * for by FOM the Dutch physics granting agency and we would like to
15  * be able to track its scientific use to convince FOM of its value
16  * for the community.
17  *
18  * This file is part of FORM.
19  *
20  * FORM is free software: you can redistribute it and/or modify it under the
21  * terms of the GNU General Public License as published by the Free Software
22  * Foundation, either version 3 of the License, or (at your option) any later
23  * version.
24  *
25  * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
26  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
27  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
28  * details.
29  *
30  * You should have received a copy of the GNU General Public License along
31  * with FORM. If not, see <http://www.gnu.org/licenses/>.
32  */
33 /* #] License : */
34 /*
35  #[ Includes : ratio.c
36 */
37 
38 #include "form3.h"
39 
40 /*
41  #] Includes :
42  #[ Ratio :
43 
44  These are the special operations regarding simple polynomials.
45  The first and most needed is the partial fractioning expansion.
46  Ratio,x1,x2,x3
47 
48  The files belonging to the ratio command serve also as a good example
49  of how to implement a new operation.
50 
51  #[ RatioFind :
52 
53  The routine that should locate the need for a ratio command.
54  If located the corresponding symbols are removed and the
55  operational parameters are loaded. A subexpression pointer
56  is inserted and the code for success is returned.
57 
58  params points at the compiler output block defined in RatioComp.
59 
60 */
61 
62 WORD RatioFind(PHEAD WORD *term, WORD *params)
63 {
64  GETBIDENTITY
65  WORD *t, *m, *r;
66  WORD x1, x2, i;
67  WORD *y1, *y2, n1 = 0, n2 = 0;
68  x1 = params[3];
69  x2 = params[4];
70  m = t = term;
71  m += *m;
72  m -= ABS(m[-1]);
73  t++;
74  if ( t < m ) do {
75  if ( *t == SYMBOL ) {
76  y1 = 0;
77  y2 = 0;
78  r = t + t[1];
79  m = t + 2;
80  do {
81  if ( *m == x1 ) { y1 = m; n1 = m[1]; }
82  else if ( *m == x2 ) { y2 = m; n2 = m[1]; }
83  m += 2;
84  } while ( m < r );
85  if ( !y1 || !y2 || ( n1 > 0 && n2 > 0 ) ) return(0);
86  m -= 2;
87  if ( y1 > y2 ) { r = y1; y1 = y2; y2 = r; }
88  *y2 = *m; y2[1] = m[1];
89  m -= 2;
90  *y1 = *m; y1[1] = m[1];
91  i = WORDDIF(m,t);
92 #if SUBEXPSIZE > 6
93 We have to revise the code for the second case.
94 #endif
95  if ( i > 2 ) { /* Subexpression fits exactly */
96  t[1] = i;
97  y1 = term+*term;
98  y2 = y1+SUBEXPSIZE-4;
99  r = m+4;
100  while ( y1 > r ) *--y2 = *--y1;
101  *m++ = SUBEXPRESSION;
102  *m++ = SUBEXPSIZE;
103  *m++ = -1;
104  *m++ = 1;
105  *m++ = DUMMYBUFFER;
106  FILLSUB(m)
107  *term += SUBEXPSIZE-4;
108  }
109  else { /* All symbols are gone. Rest has to be moved */
110  m -= 2;
111  *m++ = SUBEXPRESSION;
112  *m++ = SUBEXPSIZE;
113  *m++ = -1;
114  *m++ = 1;
115  *m++ = DUMMYBUFFER;
116  FILLSUB(m)
117  t = term;
118  t += *t;
119  *term += SUBEXPSIZE-6;
120  r = m + 6-SUBEXPSIZE;
121  do { *m++ = *r++; } while ( r < t );
122  }
123  t = AT.TMout; /* Load up the TM out array for the generator */
124  *t++ = 7;
125  *t++ = RATIO;
126  *t++ = x1;
127  *t++ = x2;
128  *t++ = params[5];
129  *t++ = n1;
130  *t++ = n2;
131  return(1);
132  }
133  t += t[1];
134  } while ( t < m );
135  return(0);
136 }
137 
138 /*
139  #] RatioFind :
140  #[ RatioGen :
141 
142  The algoritm:
143  x1^-n1*x2^n2 ==> x2 --> x1 + x3
144  x1^n1*x2^-n2 ==> x1 --> x2 - x3
145  x1^-n1*x2^-n2 ==>
146 
147  +sum(i=0,n1-1){(-1)^i*binom(n2-1+i,n2-1)
148  *x3^-(n2+i)*x1^-(n1-i)}
149  +sum(i=0,n2-1){(-1)^(n1)*binom(n1-1+i,n1-1)
150  *x3^-(n1+i)*x2^-(n2-i)}
151 
152  Actually there is an amount of arbitrariness in the first two
153  formulae and the replacement x2 -> x1 + x3 could be made 'by hand'.
154  It is better to use the nontrivial 'minimal change' formula:
155 
156  x1^-n1*x2^n2: if ( n1 >= n2 ) {
157  +sum(i=0,n2){x3^i*x1^-(n1-n2+i)*binom(n2,i)}
158  }
159  else {
160  sum(i=0,n2-n1){x2^(n2-n1-i)*x3^i*binom(n1-1+i,n1-1)}
161  +sum(i=0,n1-1){x3^(n2-i)*x1^-(n1-i)*binom(n2,i)}
162  }
163  x1^n1*x2^-n2: Same but x3 -> -x3.
164 
165  The contents of the AT.TMout/params array are:
166  length,type,x1,x2,x3,n1,n2
167 
168 */
169 
170 WORD RatioGen(PHEAD WORD *term, WORD *params, WORD num, WORD level)
171 {
172  GETBIDENTITY
173  WORD *t, *m;
174  WORD *tstops[3];
175  WORD n1, n2, i, j;
176  WORD x1,x2,x3;
177  UWORD *coef;
178  WORD ncoef, sign = 0;
179  coef = (UWORD *)AT.WorkPointer;
180  t = term;
181  tstops[2] = m = t + *t;
182  m -= ABS(m[-1]);
183  t++;
184  do {
185  if ( *t == SUBEXPRESSION && t[2] == num ) break;
186  t += t[1];
187  } while ( t < m );
188  tstops[0] = t;
189  tstops[1] = t + t[1];
190 /*
191  Copying to termout will be from term to tstop1, then the induced part
192  and finally from tstop2 to tstop3
193 
194  Now separate the various cases:
195 
196 */
197  t = params + 2;
198  x1 = *t++;
199  x2 = *t++;
200  x3 = *t++;
201  n1 = *t++;
202  n2 = *t++;
203  if ( n1 > 0 ) { /* Flip the variables and indicate -x3 */
204  n2 = -n2;
205  sign = 1;
206  i = n1; n1 = n2; n2 = i;
207  i = x1; x1 = x2; x2 = i;
208  goto PosNeg;
209  }
210  else if ( n2 > 0 ) {
211  n1 = -n1;
212 PosNeg:
213  if ( n2 <= n1 ) { /* x1 -> x2 + x3 */
214  *coef = 1;
215  ncoef = 1;
216  AT.WorkPointer = (WORD *)(coef + 1);
217  j = n2;
218  for ( i = 0; i <= n2; i++ ) {
219  if ( BinomGen(BHEAD term,level,tstops,x1,x3,n2-n1-i,i,sign&i
220  ,coef,ncoef) ) goto RatioCall;
221  if ( i < n2 ) {
222  if ( Product(coef,&ncoef,j) ) goto RatioCall;
223  if ( Quotient(coef,&ncoef,i+1) ) goto RatioCall;
224  j--;
225  AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
226  }
227  }
228  AT.WorkPointer = (WORD *)(coef);
229  return(0);
230  }
231  else {
232 /*
233  sum(i=0,n2-n1){x2^(n2-n1-i)*x3^i*binom(n1-1+i,n1-1)}
234  +sum(i=0,n1-1){x3^(n2-i)*x1^-(n1-i)*binom(n2,i)}
235 */
236  *coef = 1;
237  ncoef = 1;
238  AT.WorkPointer = (WORD *)(coef + 1);
239  j = n2 - n1;
240  for ( i = 0; i <= j; i++ ) {
241  if ( BinomGen(BHEAD term,level,tstops,x2,x3,n2-n1-i,i,sign&i
242  ,coef,ncoef) ) goto RatioCall;
243  if ( i < j ) {
244  if ( Product(coef,&ncoef,n1+i) ) goto RatioCall;
245  if ( Quotient(coef,&ncoef,i+1) ) goto RatioCall;
246  AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
247  }
248  }
249  *coef = 1;
250  ncoef = 1;
251  AT.WorkPointer = (WORD *)(coef + 1);
252  j = n1-1;
253  for ( i = 0; i <= j; i++ ) {
254  if ( BinomGen(BHEAD term,level,tstops,x1,x3,i-n1,n2-i,sign&(n2-i)
255  ,coef,ncoef) ) goto RatioCall;
256  if ( i < j ) {
257  if ( Product(coef,&ncoef,n2-i) ) goto RatioCall;
258  if ( Quotient(coef,&ncoef,i+1) ) goto RatioCall;
259  AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
260  }
261  }
262  AT.WorkPointer = (WORD *)(coef);
263  return(0);
264  }
265  }
266  else {
267  n2 = -n2;
268  n1 = -n1;
269 /*
270  +sum(i=0,n1-1){(-1)^i*binom(n2-1+i,n2-1)
271  *x3^-(n2+i)*x1^-(n1-i)}
272  +sum(i=0,n2-1){(-1)^(n1)*binom(n1-1+i,n1-1)
273  *x3^-(n1+i)*x2^-(n2-i)}
274 */
275  *coef = 1;
276  ncoef = 1;
277  AT.WorkPointer = (WORD *)(coef + 1);
278  j = n1-1;
279  for ( i = 0; i <= j; i++ ) {
280  if ( BinomGen(BHEAD term,level,tstops,x1,x3,i-n1,-n2-i,i&1
281  ,coef,ncoef) ) goto RatioCall;
282  if ( i < j ) {
283  if ( Product(coef,&ncoef,n2+i) ) goto RatioCall;
284  if ( Quotient(coef,&ncoef,i+1) ) goto RatioCall;
285  AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
286  }
287  }
288  *coef = 1;
289  ncoef = 1;
290  AT.WorkPointer = (WORD *)(coef + 1);
291  j = n2-1;
292  for ( i = 0; i <= j; i++ ) {
293  if ( BinomGen(BHEAD term,level,tstops,x2,x3,i-n2,-n1-i,n1&1
294  ,coef,ncoef) ) goto RatioCall;
295  if ( i < j ) {
296  if ( Product(coef,&ncoef,n1+i) ) goto RatioCall;
297  if ( Quotient(coef,&ncoef,i+1) ) goto RatioCall;
298  AT.WorkPointer = (WORD *)(coef + ABS(ncoef));
299  }
300  }
301  AT.WorkPointer = (WORD *)(coef);
302  return(0);
303  }
304 
305 RatioCall:
306  MLOCK(ErrorMessageLock);
307  MesCall("RatioGen");
308  MUNLOCK(ErrorMessageLock);
309  SETERROR(-1)
310 }
311 
312 /*
313  #] RatioGen :
314  #[ BinomGen :
315 
316  Routine for the generation of terms in a binomialtype expansion.
317 
318 */
319 
320 WORD BinomGen(PHEAD WORD *term, WORD level, WORD **tstops, WORD x1, WORD x2,
321  WORD pow1, WORD pow2, WORD sign, UWORD *coef, WORD ncoef)
322 {
323  GETBIDENTITY
324  WORD *t, *r;
325  WORD *termout;
326  WORD k;
327  termout = AT.WorkPointer;
328  t = termout;
329  r = term;
330  do { *t++ = *r++; } while ( r < tstops[0] );
331  *t++ = SYMBOL;
332  if ( pow2 == 0 ) {
333  if ( pow1 == 0 ) t--;
334  else { *t++ = 4; *t++ = x1; *t++ = pow1; }
335  }
336  else if ( pow1 == 0 ) {
337  *t++ = 4; *t++ = x2; *t++ = pow2;
338  }
339  else {
340  *t++ = 6; *t++ = x1; *t++ = pow1; *t++ = x2; *t++ = pow2;
341  }
342  *t++ = LNUMBER;
343  *t++ = ABS(ncoef) + 3;
344  *t = ncoef;
345  if ( sign ) *t = -*t;
346  t++;
347  ncoef = ABS(ncoef);
348  for ( k = 0; k < ncoef; k++ ) *t++ = coef[k];
349  r = tstops[1];
350  do { *t++ = *r++; } while ( r < tstops[2] );
351  *termout = WORDDIF(t,termout);
352  AT.WorkPointer = t;
353  if ( AT.WorkPointer > AT.WorkTop ) {
354  MLOCK(ErrorMessageLock);
355  MesWork();
356  MUNLOCK(ErrorMessageLock);
357  return(-1);
358  }
359  *AN.RepPoint = 1;
360  AR.expchanged = 1;
361  if ( Generator(BHEAD termout,level) ) {
362  MLOCK(ErrorMessageLock);
363  MesCall("BinomGen");
364  MUNLOCK(ErrorMessageLock);
365  SETERROR(-1)
366  }
367  AT.WorkPointer = termout;
368  return(0);
369 }
370 
371 /*
372  #] BinomGen :
373  #] Ratio :
374  #[ Sum :
375  #[ DoSumF1 :
376 
377  Routine expands a sum_ function.
378  Its arguments are:
379  The term in which the function occurs.
380  The parameter list:
381  length of parameter field
382  function number (SUMNUM1)
383  number of the symbol that is loop parameter
384  min value
385  max value
386  increment
387  the number of the subexpression to be removed
388  the level in the generation tree.
389 
390  Note that the insertion of the loop parameter in the argument
391  is done via the regular wildcard substitution mechanism.
392 
393 */
394 
395 WORD DoSumF1(PHEAD WORD *term, WORD *params, WORD replac, WORD level)
396 {
397  GETBIDENTITY
398  WORD *termout, *t, extractbuff = AT.TMbuff;
399  WORD isum, ival, iinc;
400  LONG from;
401  CBUF *C;
402  ival = params[3];
403  iinc = params[5];
404  if ( ( iinc > 0 && params[4] >= ival )
405  || ( iinc < 0 && params[4] <= ival ) ) {
406  isum = (params[4] - ival)/iinc + 1;
407  }
408  else return(0);
409  termout = AT.WorkPointer;
410  AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
411  if ( AT.WorkPointer > AT.WorkTop ) {
412  MLOCK(ErrorMessageLock);
413  MesWork();
414  MUNLOCK(ErrorMessageLock);
415  return(-1);
416  }
417  t = term + 1;
418  while ( *t != SUBEXPRESSION || t[2] != replac || t[4] != extractbuff )
419  t += t[1];
420  C = cbuf+t[4];
421  t += SUBEXPSIZE;
422  if ( params[2] < 0 ) {
423  while ( *t != INDTOIND || t[2] != -params[2] ) t += t[1];
424  *t = INDTOIND;
425  }
426  else {
427  while ( *t > SYMTOSUB || t[2] != params[2] ) t += t[1];
428  *t = SYMTONUM;
429  }
430  do {
431  t[3] = ival;
432  from = C->rhs[replac] - C->Buffer;
433  while ( C->Buffer[from] ) {
434  if ( InsertTerm(BHEAD term,replac,extractbuff,C->Buffer+from,termout,0) < 0 ) goto SumF1Call;
435  AT.WorkPointer = termout + *termout;
436  if ( Generator(BHEAD termout,level) < 0 ) goto SumF1Call;
437  from += C->Buffer[from];
438  }
439  ival += iinc;
440  } while ( --isum > 0 );
441  AT.WorkPointer = termout;
442  return(0);
443 SumF1Call:
444  MLOCK(ErrorMessageLock);
445  MesCall("DoSumF1");
446  MUNLOCK(ErrorMessageLock);
447  SETERROR(-1)
448 }
449 
450 /*
451  #] DoSumF1 :
452  #[ Glue :
453 
454  Routine multiplies two terms. The second term is subject
455  to the wildcard substitutions in sub.
456  Output in the first term. This routine is a variation on
457  the routine InsertTerm.
458 
459 */
460 
461 WORD Glue(PHEAD WORD *term1, WORD *term2, WORD *sub, WORD insert)
462 {
463  GETBIDENTITY
464  UWORD *coef;
465  WORD ncoef, *t, *t1, *t2, i, nc2, nc3, old, newer;
466  coef = (UWORD *)(TermMalloc("Glue"));
467  t = term1;
468  t += *t;
469  i = t[-1];
470  t -= ABS(i);
471  old = WORDDIF(t,term1);
472  ncoef = REDLENG(i);
473  if ( i < 0 ) i = -i;
474  i--;
475  t1 = t;
476  t2 = (WORD *)coef;
477  while ( --i >= 0 ) *t2++ = *t1++;
478  i = *--t;
479  nc2 = WildFill(BHEAD t,term2,sub);
480  *t = i;
481  t += nc2;
482  nc2 = t[-1];
483  t -= ABS(nc2);
484  newer = WORDDIF(t,term1);
485  if ( MulRat(BHEAD (UWORD *)t,REDLENG(nc2),coef,ncoef,(UWORD *)t,&nc3) ) {
486  MLOCK(ErrorMessageLock);
487  MesCall("Glue");
488  MUNLOCK(ErrorMessageLock);
489  TermFree(coef,"Glue");
490  SETERROR(-1)
491  }
492  i = (ABS(nc3))<<1;
493  t += i++;
494  *t++ = (nc3 >= 0)?i:-i;
495  *term1 = WORDDIF(t,term1);
496 /*
497  Switch the new piece with the old tail, so that noncommuting
498  variables get into their proper spot.
499 */
500  i = old - insert;
501  t1 = t;
502  t2 = term1+insert;
503  NCOPY(t1,t2,i);
504  i = newer - old;
505  t1 = term1+insert;
506  t2 = term1+old;
507  NCOPY(t1,t2,i);
508  t2 = t;
509  i = old - insert;
510  NCOPY(t1,t2,i);
511  TermFree(coef,"Glue");
512  return(0);
513 }
514 
515 /*
516  #] Glue :
517  #[ DoSumF2 :
518 */
519 
520 WORD DoSumF2(PHEAD WORD *term, WORD *params, WORD replac, WORD level)
521 {
522  GETBIDENTITY
523  WORD *termout, *t, *from, *sub, *to, extractbuff = AT.TMbuff;
524  WORD isum, ival, iinc, insert, i;
525  CBUF *C;
526  ival = params[3];
527  iinc = params[5];
528  if ( ( iinc > 0 && params[4] >= ival )
529  || ( iinc < 0 && params[4] <= ival ) ) {
530  isum = (params[4] - ival)/iinc + 1;
531  }
532  else return(0);
533  termout = AT.WorkPointer;
534  AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
535  if ( AT.WorkPointer > AT.WorkTop ) {
536  MLOCK(ErrorMessageLock);
537  MesWork();
538  MUNLOCK(ErrorMessageLock);
539  return(-1);
540  }
541  t = term + 1;
542  while ( *t != SUBEXPRESSION || t[2] != replac || t[4] != extractbuff ) t += t[1];
543  insert = WORDDIF(t,term);
544 
545  from = term;
546  to = termout;
547  while ( from < t ) *to++ = *from++;
548  from += t[1];
549  sub = term + *term;
550  while ( from < sub ) *to++ = *from++;
551  *termout -= t[1];
552 
553  sub = t;
554  C = cbuf+t[4];
555  t += SUBEXPSIZE;
556  if ( params[2] < 0 ) {
557  while ( *t != INDTOIND || t[2] != -params[2] ) t += t[1];
558  *t = INDTOIND;
559  }
560  else {
561  while ( *t > SYMTOSUB || t[2] != params[2] ) t += t[1];
562  *t = SYMTONUM;
563  }
564  t[3] = ival;
565  for(;;) {
566  AT.WorkPointer = termout + *termout;
567  to = AT.WorkPointer;
568  if ( ( to + *termout ) > AT.WorkTop ) {
569  MLOCK(ErrorMessageLock);
570  MesWork();
571  MUNLOCK(ErrorMessageLock);
572  return(-1);
573  }
574  from = termout;
575  i = *termout;
576  NCOPY(to,from,i);
577  from = AT.WorkPointer;
578  AT.WorkPointer = to;
579  if ( Generator(BHEAD from,level) < 0 ) goto SumF2Call;
580  if ( --isum <= 0 ) break;
581  ival += iinc;
582  t[3] = ival;
583  if ( Glue(BHEAD termout,C->rhs[replac],sub,insert) < 0 ) goto SumF2Call;
584  }
585  AT.WorkPointer = termout;
586  return(0);
587 SumF2Call:
588  MLOCK(ErrorMessageLock);
589  MesCall("DoSumF2");
590  MUNLOCK(ErrorMessageLock);
591  SETERROR(-1)
592 }
593 
594 /*
595  #] DoSumF2 :
596  #] Sum :
597  #[ GCDfunction :
598  #[ GCDfunction :
599 */
600 
601 typedef struct {
602  WORD *buffer;
603  DOLLARS dollar;
604  LONG size;
605  int type;
606  int dummy;
607 } ARGBUFFER;
608 
609 int GCDfunction(PHEAD WORD *term,WORD level)
610 {
611  GETBIDENTITY
612  WORD *t, *tstop, *tf, *termout, *tin, *tout, *m, *mnext, *mstop, *mm;
613  int todo, i, ii, j, istart, sign = 1, action = 0;
614  WORD firstshort = 0, firstvalue = 0, gcdisone = 0, mlength, tlength, newlength;
615  WORD totargs = 0, numargs, *mh, oldval1, *g, *gcdout = 0;
616  WORD *arg1, *arg2;
617  UWORD x1,x2,x3;
618  LONG args;
619 #if ( FUNHEAD > 4 )
620  WORD sh[FUNHEAD+5];
621 #else
622  WORD sh[9];
623 #endif
624  DOLLARS d;
625  ARGBUFFER *abuf = 0, ab;
626 /*
627  #[ Find Function. Count arguments :
628 
629  First find the proper function
630 */
631  t = term + *term; tlength = t[-1];
632  tstop = t - ABS(tlength);
633  t = term + 1;
634  while ( t < tstop ) {
635  if ( *t != GCDFUNCTION ) { t += t[1]; continue; }
636  todo = 1; totargs = 0;
637  tf = t + FUNHEAD;
638  while ( tf < t + t[1] ) {
639  totargs++;
640  if ( *tf > 0 && tf[1] != 0 ) todo = 0;
641  NEXTARG(tf);
642  }
643  if ( todo ) break;
644  t += t[1];
645  }
646  if ( t >= tstop ) {
647  MLOCK(ErrorMessageLock);
648  MesPrint("Internal error. Indicated gcd_ function not encountered.");
649  MUNLOCK(ErrorMessageLock);
650  Terminate(-1);
651  }
652  WantAddPointers(totargs);
653  args = AT.pWorkPointer; AT.pWorkPointer += totargs;
654 /*
655  #] Find Function. Count arguments :
656  #[ Do short arguments :
657 
658  The function we need, in agreement with TestSub, is now in t
659  Make first a compilation of the short arguments (except $-s and expressions)
660  to see whether we need to do much work.
661  This means that after this scan we can ignore all short arguments with
662  the exception of unevaluated $-s and expressions.
663 */
664  numargs = 0;
665  firstshort = 0;
666  tf = t + FUNHEAD;
667  while ( tf < t + t[1] ) {
668  if ( *tf == -SNUMBER && tf[1] == 0 ) { NEXTARG(tf); continue; }
669  if ( *tf > 0 || *tf == -DOLLAREXPRESSION || *tf == -EXPRESSION ) {
670  AT.pWorkSpace[args+numargs++] = tf;
671  NEXTARG(tf); continue;
672  }
673  if ( firstshort == 0 ) {
674  firstshort = *tf;
675  if ( *tf <= -FUNCTION ) { firstvalue = -(*tf); }
676  else { firstvalue = tf[1]; }
677  NEXTARG(tf);
678  continue;
679  }
680  else if ( *tf != firstshort ) {
681  if ( *tf != -INDEX && *tf != -VECTOR && *t != -MINVECTOR ) {
682  gcdisone = 1; break;
683  }
684  if ( firstshort != -INDEX && firstshort != -VECTOR && firstshort != -MINVECTOR ) {
685  gcdisone = 1; break;
686  }
687  if ( tf[1] != firstvalue ) {
688  gcdisone = 1; break;
689  }
690  if ( *t == -MINVECTOR ) { firstshort = -VECTOR; }
691  if ( firstshort == -MINVECTOR ) { firstshort = -VECTOR; }
692  }
693  else if ( *tf > -FUNCTION && *tf != -SNUMBER && tf[1] != firstvalue ) {
694  gcdisone = 1;
695  break;
696  }
697  if ( *tf == -SNUMBER && firstvalue != tf[1] ) {
698 /*
699  make a new firstvalue which is gcd_(firstvalue,tf[1])
700 */
701  if ( firstvalue == 1 || tf[1] == 1 ) { gcdisone = 1; break; }
702  if ( firstvalue < 0 && tf[1] < 0 ) {
703  x1 = -firstvalue; x2 = -tf[1]; sign = -1;
704  }
705  else {
706  x1 = ABS(firstvalue); x2 = ABS(tf[1]); sign = 1;
707  }
708  while ( ( x3 = x1%x2 ) != 0 ) { x1 = x2; x2 = x3; }
709  firstvalue = ((WORD)x2)*sign;
710  if ( firstvalue == 1 ) { gcdisone = 1; break; }
711  }
712  NEXTARG(tf);
713  }
714  termout = AT.WorkPointer;
715  AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
716  if ( AT.WorkPointer > AT.WorkTop ) {
717  MLOCK(ErrorMessageLock);
718  MesWork();
719  MUNLOCK(ErrorMessageLock);
720  return(-1);
721  }
722 /*
723  #] Do short arguments :
724  #[ Do trivial GCD :
725 
726  Copy head
727 */
728  i = t - term; tin = term; tout = termout;
729  NCOPY(tout,tin,i);
730  if ( gcdisone || ( firstshort == -SNUMBER && firstvalue == 1 ) ) {
731  sign = 1;
732 gcdone:
733  tin += t[1]; tstop = term + *term;
734  while ( tin < tstop ) *tout++ = *tin++;
735  *termout = tout - termout;
736  if ( sign < 0 ) tout[-1] = -tout[-1];
737  AT.WorkPointer = tout;
738  if ( Generator(BHEAD termout,level) < 0 ) goto CalledFrom;
739  AT.WorkPointer = termout;
740  AT.pWorkPointer = args;
741  return(0);
742  }
743 /*
744  #] Do trivial GCD :
745  #[ Do short argument GCD :
746 */
747  if ( numargs == 0 ) { /* basically we are done */
748 doshort:
749  sign = 1;
750  if ( firstshort == 0 ) goto gcdone;
751  if ( firstshort == -SNUMBER ) {
752  *tout++ = SNUMBER; *tout++ = 4; *tout++ = firstvalue; *tout++ = 1;
753  goto gcdone;
754  }
755  else if ( firstshort == -SYMBOL ) {
756  *tout++ = SYMBOL; *tout++ = 4; *tout++ = firstvalue; *tout++ = 1;
757  goto gcdone;
758  }
759  else if ( firstshort == -VECTOR || firstshort == -INDEX ) {
760  *tout++ = INDEX; *tout++ = 3; *tout++ = firstvalue; goto gcdone;
761  }
762  else if ( firstshort == -MINVECTOR ) {
763  sign = -1;
764  *tout++ = INDEX; *tout++ = 3; *tout++ = firstvalue; goto gcdone;
765  }
766  else if ( firstshort <= -FUNCTION ) {
767  *tout++ = firstvalue; *tout++ = FUNHEAD; FILLFUN(tout);
768  goto gcdone;
769  }
770  else {
771  MLOCK(ErrorMessageLock);
772  MesPrint("Internal error. Illegal short argument in GCDfunction.");
773  MUNLOCK(ErrorMessageLock);
774  Terminate(-1);
775  }
776  }
777 /*
778  #] Do short argument GCD :
779  #[ Convert short argument :
780 
781  Now we allocate space for the arguments in general notation.
782  First the special one if there were short arguments
783 */
784  if ( firstshort ) {
785  switch ( firstshort ) {
786  case -SNUMBER:
787  sh[0] = 4; sh[1] = ABS(firstvalue); sh[2] = 1;
788  if ( firstvalue < 0 ) sh[3] = -3;
789  else sh[3] = 3;
790  sh[4] = 0;
791  break;
792  case -MINVECTOR:
793  case -VECTOR:
794  case -INDEX:
795  sh[0] = 8; sh[1] = INDEX; sh[2] = 3; sh[3] = firstvalue;
796  sh[4] = 1; sh[5] = 1;
797  if ( firstshort == -MINVECTOR ) sh[6] = -3;
798  else sh[6] = 3;
799  sh[7] = 0;
800  break;
801  case -SYMBOL:
802  sh[0] = 8; sh[1] = SYMBOL; sh[2] = 4; sh[3] = firstvalue; sh[4] = 1;
803  sh[5] = 1; sh[6] = 1; sh[7] = 3; sh[8] = 0;
804  break;
805  default:
806  sh[0] = FUNHEAD+4; sh[1] = firstshort; sh[2] = FUNHEAD;
807  for ( i = 2; i < FUNHEAD; i++ ) sh[i+1] = 0;
808  sh[FUNHEAD+1] = 1; sh[FUNHEAD+2] = 1; sh[FUNHEAD+3] = 3; sh[FUNHEAD+4] = 0;
809  break;
810  }
811  }
812 /*
813  #] Convert short argument :
814  #[ Sort arguments :
815 
816  Now we should sort the arguments in a way that the dollars and the
817  expressions come last. That way we may never need them.
818 */
819  for ( i = 1; i < numargs; i++ ) {
820  for ( ii = i; ii > 0; ii-- ) {
821  arg1 = AT.pWorkSpace[args+ii];
822  arg2 = AT.pWorkSpace[args+ii-1];
823  if ( *arg1 < 0 ) {
824  if ( *arg2 < 0 ) {
825  if ( *arg1 == -EXPRESSION ) break;
826  if ( *arg2 == -DOLLAREXPRESSION ) break;
827  AT.pWorkSpace[args+ii] = arg2;
828  AT.pWorkSpace[args+ii-1] = arg1;
829  }
830  else break;
831  }
832  else if ( *arg2 < 0 ) {
833  AT.pWorkSpace[args+ii] = arg2;
834  AT.pWorkSpace[args+ii-1] = arg1;
835  }
836  else {
837  if ( *arg1 > *arg2 ) {
838  AT.pWorkSpace[args+ii] = arg2;
839  AT.pWorkSpace[args+ii-1] = arg1;
840  }
841  else break;
842  }
843  }
844  }
845 /*
846  #] Sort arguments :
847  #[ There is a single term argument :
848 */
849  if ( firstshort ) {
850  mh = sh; istart = 0;
851 oneterm:;
852  for ( i = istart; i < numargs; i++ ) {
853  arg1 = AT.pWorkSpace[args+i];
854  if ( *arg1 > 0 ) {
855  oldval1 = arg1[*arg1]; arg1[*arg1] = 0;
856  m = arg1+ARGHEAD;
857  while ( *m ) {
858  GCDterms(BHEAD mh,m,mh); m += *m;
859  if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
860  gcdisone = 1; sign = 1; arg1[*arg1] = oldval1; goto gcdone;
861  }
862  }
863  arg1[*arg1] = oldval1;
864  }
865  else if ( *arg1 == -DOLLAREXPRESSION ) {
866  if ( ( d = DolToTerms(BHEAD arg1[1]) ) != 0 ) {
867  m = d->where;
868  while ( *m ) {
869  GCDterms(BHEAD mh,m,mh); m += *m;
870  if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
871  gcdisone = 1; sign = 1;
872  if ( d->factors ) M_free(d->factors,"Dollar factors");
873  M_free(d,"Copy of dollar variable"); goto gcdone;
874  }
875  }
876  if ( d->factors ) M_free(d->factors,"Dollar factors");
877  M_free(d,"Copy of dollar variable");
878  }
879  }
880  else {
881  mm = CreateExpression(BHEAD arg1[1]);
882  m = mm;
883  while ( *m ) {
884  GCDterms(BHEAD mh,m,mh); m += *m;
885  if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
886  gcdisone = 1; sign = 1; M_free(mm,"CreateExpression"); goto gcdone;
887  }
888  }
889  M_free(mm,"CreateExpression");
890  }
891  }
892  if ( firstshort ) {
893  if ( mh[0] == 4 ) {
894  firstshort = -SNUMBER; firstvalue = mh[1] * (mh[3]/3);
895  }
896  else if ( mh[1] == SYMBOL ) {
897  firstshort = -SYMBOL; firstvalue = mh[3];
898  }
899  else if ( mh[1] == INDEX ) {
900  firstshort = -INDEX; firstvalue = mh[3];
901  if ( mh[6] == -3 ) firstshort = -MINVECTOR;
902  }
903  else if ( mh[1] >= FUNCTION ) {
904  firstshort = -mh[1]; firstvalue = mh[1];
905  }
906  goto doshort;
907  }
908  else {
909 /*
910  We have a GCD that is only a single term.
911  Paste it in and combine the coefficients.
912 */
913  mh[mh[0]] = 0;
914  mm = mh;
915  goto multiterms;
916  }
917  }
918 /*
919  Now we have only regular arguments.
920  But some have not yet been expanded.
921  Check whether there are proper long arguments and if so if there is
922  one with just a single term
923 */
924  for ( i = 0; i < numargs; i++ ) {
925  arg1 = AT.pWorkSpace[args+i];
926  if ( *arg1 > 0 && arg1[ARGHEAD]+ARGHEAD == *arg1 ) {
927 /*
928  We have an argument with a single term
929 */
930  if ( i != 0 ) {
931  arg2 = AT.pWorkSpace[args];
932  AT.pWorkSpace[args] = arg1;
933  AT.pWorkSpace[args+1] = arg2;
934  }
935  m = mh = AT.WorkPointer;
936  mm = arg1+ARGHEAD; i = *mm;
937  NCOPY(m,mm,i);
938  AT.WorkPointer = m;
939  istart = 1;
940  goto oneterm;
941  }
942  }
943 /*
944  #] There is a single term argument :
945  #[ Expand $ and expr :
946 
947  We have: 1: regular multiterm arguments
948  2: dollars
949  3: expressions.
950  The sum of them is numargs. Their addresses are in args. The problem is
951  that expansion will lead to allocations that we have to return and all
952  these allocations are different in nature.
953 */
954  action = 1;
955  abuf = (ARGBUFFER *)Malloc1(numargs*sizeof(ARGBUFFER),"argbuffer");
956  for ( i = 0; i < numargs; i++ ) {
957  arg1 = AT.pWorkSpace[args+i];
958  if ( *arg1 > 0 ) {
959  m = (WORD *)Malloc1(*arg1*sizeof(WORD),"argbuffer type 0");
960  abuf[i].buffer = m;
961  abuf[i].type = 0;
962  mm = arg1+ARGHEAD;
963  j = *arg1-ARGHEAD;
964  abuf[i].size = j;
965  NCOPY(m,mm,j);
966  *m = 0;
967  }
968  else if ( *arg1 == -DOLLAREXPRESSION ) {
969  d = DolToTerms(BHEAD arg1[1]);
970  abuf[i].buffer = d->where;
971  abuf[i].type = 1;
972  abuf[i].dollar = d;
973  m = abuf[i].buffer; while ( *m ) m+= *m;
974  abuf[i].size = m-abuf[i].buffer;
975  }
976  else if ( *arg1 == -EXPRESSION ) {
977  abuf[i].buffer = CreateExpression(BHEAD arg1[1]);
978  abuf[i].type = 2;
979  m = abuf[i].buffer; while ( *m ) m+= *m;
980  abuf[i].size = m-abuf[i].buffer;
981  }
982  else {
983  MLOCK(ErrorMessageLock);
984  MesPrint("What argument is this?");
985  MUNLOCK(ErrorMessageLock);
986  goto CalledFrom;
987  }
988  }
989  for ( i = 0; i < numargs; i++ ) {
990  arg1 = abuf[i].buffer;
991  if ( arg1[*arg1] == 0 ) {
992 /*
993  After expansion there is an argument with a single term
994 */
995  ab = abuf[i]; abuf[i] = abuf[0]; abuf[0] = ab;
996  mh = abuf[0].buffer;
997  for ( j = 1; j < numargs; j++ ) {
998  m = abuf[j].buffer;
999  while ( *m ) {
1000  GCDterms(BHEAD mh,m,mh); m += *m;
1001  if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
1002  gcdisone = 1; sign = 1; break;
1003  }
1004  }
1005  if ( *m ) break;
1006  }
1007  mm = mh + *mh; if ( mm[-1] < 0 ) { sign = -1; mm[-1] = -mm[-1]; }
1008  mstop = mm - mm[-1]; m = mh+1; mlength = mm[-1];
1009  while ( tin < t ) *tout++ = *tin++;
1010  while ( m < mstop ) *tout++ = *m++;
1011  tin += tin[1];
1012  while ( tin < tstop ) *tout++ = *tin++;
1013  tlength = REDLENG(tlength);
1014  mlength = REDLENG(mlength);
1015  if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)mstop,mlength,
1016  (UWORD *)tout,&newlength) < 0 ) goto CalledFrom;
1017  mlength = INCLENG(newlength);
1018  tout += ABS(mlength);
1019  tout[-1] = mlength*sign;
1020  *termout = tout - termout;
1021  AT.WorkPointer = tout;
1022  if ( Generator(BHEAD termout,level) < 0 ) goto CalledFrom;
1023  goto cleanup;
1024  }
1025  }
1026 /*
1027  There are only arguments with more than one term.
1028  We order them by size to make the computations as easy as possible.
1029 */
1030  for ( i = 1; i < numargs; i++ ) {
1031  for ( ii = i; ii > 0; ii-- ) {
1032  if ( abuf[ii-1].size <= abuf[ii].size ) break;
1033  ab = abuf[ii-1]; abuf[ii-1] = abuf[ii]; abuf[ii] = ab;
1034  }
1035  }
1036 /*
1037  #] Expand $ and expr :
1038  #[ Multiterm subexpressions :
1039 */
1040  gcdout = abuf[0].buffer;
1041  for ( i = 1; i < numargs; i++ ) {
1042  g = GCDfunction3(BHEAD gcdout,abuf[i].buffer);
1043  if ( gcdout != abuf[0].buffer ) M_free(gcdout,"gcdout");
1044  gcdout = g;
1045  if ( gcdout[*gcdout] == 0 && gcdout[0] == 4 && gcdout[1] == 1
1046  && gcdout[2] == 1 && gcdout[3] == 3 ) break;
1047  }
1048  mm = gcdout;
1049 multiterms:;
1050  tlength = REDLENG(tlength);
1051  while ( *mm ) {
1052  tin = term; tout = termout; while ( tin < t ) *tout++ = *tin++;
1053  tin += t[1];
1054  mnext = mm + *mm; mlength = mnext[-1]; mstop = mnext - ABS(mlength);
1055  mm++;
1056  while ( mm < mstop ) *tout++ = *mm++;
1057  while ( tin < tstop ) *tout++ = *tin++;
1058  mlength = REDLENG(mlength);
1059  if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)mm,mlength,
1060  (UWORD *)tout,&newlength) < 0 ) goto CalledFrom;
1061  mlength = INCLENG(newlength);
1062  tout += ABS(mlength);
1063  tout[-1] = mlength;
1064  *termout = tout - termout;
1065  AT.WorkPointer = tout;
1066  if ( Generator(BHEAD termout,level) < 0 ) goto CalledFrom;
1067  mm = mnext; /* next term */
1068  }
1069  if ( action && ( gcdout != abuf[0].buffer ) ) M_free(gcdout,"gcdout");
1070 /*
1071  #] Multiterm subexpressions :
1072  #[ Cleanup :
1073 */
1074 cleanup:;
1075  if ( action ) {
1076  for ( i = 0; i < numargs; i++ ) {
1077  if ( abuf[i].type == 0 ) { M_free(abuf[i].buffer,"argbuffer type 0"); }
1078  else if ( abuf[i].type == 1 ) {
1079  d = abuf[i].dollar;
1080  if ( d->factors ) M_free(d->factors,"Dollar factors");
1081  M_free(d,"Copy of dollar variable");
1082  }
1083  else if ( abuf[i].type == 2 ) { M_free(abuf[i].buffer,"CreateExpression"); }
1084  }
1085  M_free(abuf,"argbuffer");
1086  }
1087 /*
1088  #] Cleanup :
1089 */
1090  AT.pWorkPointer = args;
1091  AT.WorkPointer = termout;
1092  return(0);
1093 
1094 CalledFrom:
1095  MLOCK(ErrorMessageLock);
1096  MesCall("GCDfunction");
1097  MUNLOCK(ErrorMessageLock);
1098  SETERROR(-1)
1099  return(-1);
1100 }
1101 
1102 /*
1103  #] GCDfunction :
1104  #[ GCDfunction3 :
1105 
1106  Finds the GCD of the two arguments which are buffers with terms.
1107  In principle the first buffer can have only one term.
1108 
1109  If both buffers have more than one term, we need to replace all
1110  non-symbolic objects by generated symbols and substitute that back
1111  afterwards. The rest we leave to the powerful routines.
1112  Philosophical problem: What do we do with GCD_(x/z+y,x+y*z) ?
1113 
1114  Method:
1115  If we have only negative powers of z and no positive powers we let
1116  the EXTRASYMBOLS do their job. When mixed, multiply the arguments with
1117  the negative powers with enough powers of z to eliminate the negative powers.
1118  The DENOMINATOR function is always eliminated with the mechanism as we
1119  cannot tell whether there are positive powers of its contents.
1120 */
1121 
1122 WORD *GCDfunction3(PHEAD WORD *in1, WORD *in2)
1123 {
1124  GETBIDENTITY
1125  WORD oldsorttype = AR.SortType;
1126  WORD *t, *tt, *gcdout, *term1, *term2, *confree1, *confree2, *gcdout1, *proper1, *proper2;
1127  int i, actionflag1, actionflag2;
1128  WORD startebuf = cbuf[AT.ebufnum].numrhs;
1129  if ( in2[*in2] == 0 ) { t = in1; in1 = in2; in2 = t; }
1130  if ( in1[*in1] == 0 ) { /* First input with only one term */
1131  gcdout = (WORD *)Malloc1((*in1+1)*sizeof(WORD),"gcdout");
1132  i = *in1; t = gcdout; tt = in1; NCOPY(t,tt,i); *t = 0;
1133  t = in2;
1134  while ( *t ) {
1135  GCDterms(BHEAD gcdout,t,gcdout);
1136  if ( gcdout[0] == 4 && gcdout[1] == 1
1137  && gcdout[2] == 1 && gcdout[3] == 3 ) break;
1138  t += *t;
1139  }
1140  return(gcdout);
1141  }
1142 /*
1143  We need to take out the content from the two expressions
1144  and determine their GCD. This plays with the negative powers!
1145 */
1146  AR.SortType = SORTHIGHFIRST;
1147  term1 = TermMalloc("GCDfunction3-a");
1148  term2 = TermMalloc("GCDfunction3-b");
1149  confree1 = TakeContent(BHEAD in1,term1);
1150  confree2 = TakeContent(BHEAD in2,term2);
1151  GCDterms(BHEAD term1,term2,term1);
1152  TermFree(term2,"GCDfunction3-b");
1153 /*
1154  Now we have to replace all non-symbols and symbols to a negative power
1155  by extra symbols.
1156 */
1157  if ( ( proper1 = PutExtraSymbols(BHEAD confree1,startebuf,&actionflag1) ) == 0 ) goto CalledFrom;
1158  if ( confree1 != in1 ) M_free(confree1,"TakeContent");
1159  if ( ( proper2 = PutExtraSymbols(BHEAD confree2,startebuf,&actionflag2) ) == 0 ) goto CalledFrom;
1160  if ( confree2 != in2 ) M_free(confree2,"TakeContent");
1161 /*
1162  And now the real work:
1163 */
1164  gcdout1 = poly_gcd(BHEAD proper1,proper2);
1165  M_free(proper1,"PutExtraSymbols");
1166  M_free(proper2,"PutExtraSymbols");
1167 
1168  AR.SortType = oldsorttype;
1169  if ( actionflag1 || actionflag2 ) {
1170  if ( ( gcdout = TakeExtraSymbols(BHEAD gcdout1,startebuf) ) == 0 ) goto CalledFrom;
1171  M_free(gcdout1,"gcdout");
1172  }
1173  else {
1174  gcdout = gcdout1;
1175  }
1176 
1177  cbuf[AT.ebufnum].numrhs = startebuf;
1178 /*
1179  Now multiply gcdout by term1
1180 */
1181  if ( term1[0] != 4 || term1[3] != 3 || term1[1] != 1 || term1[2] != 1 ) {
1182  if ( ( gcdout1 = MultiplyWithTerm(BHEAD gcdout,term1) ) == 0 ) goto CalledFrom;
1183  M_free(gcdout,"gcdout");
1184  gcdout = gcdout1;
1185  }
1186  TermFree(term1,"GCDfunction3-a");
1187  return(gcdout);
1188 CalledFrom:
1189  MLOCK(ErrorMessageLock);
1190  MesCall("GCDfunction3");
1191  MUNLOCK(ErrorMessageLock);
1192  return(0);
1193 }
1194 
1195 /*
1196  #] GCDfunction3 :
1197  #[ PutExtraSymbols :
1198 */
1199 
1200 WORD *PutExtraSymbols(PHEAD WORD *in,WORD startebuf,int *actionflag)
1201 {
1202  WORD *termout = AT.WorkPointer;
1203  int action;
1204  *actionflag = 0;
1205  NewSort(BHEAD0);
1206  while ( *in ) {
1207  if ( ( action = LocalConvertToPoly(BHEAD in,termout,startebuf,0) ) < 0 ) {
1208  LowerSortLevel();
1209  goto CalledFrom;
1210  }
1211  if ( action > 0 ) *actionflag = 1;
1212  StoreTerm(BHEAD termout);
1213  in += *in;
1214  }
1215  if ( EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 ) goto CalledFrom;
1216  return(termout);
1217 CalledFrom:
1218  MLOCK(ErrorMessageLock);
1219  MesCall("PutExtraSymbols");
1220  MUNLOCK(ErrorMessageLock);
1221  return(0);
1222 }
1223 
1224 /*
1225  #] PutExtraSymbols :
1226  #[ TakeExtraSymbols :
1227 */
1228 
1229 WORD *TakeExtraSymbols(PHEAD WORD *in,WORD startebuf)
1230 {
1231  CBUF *C = cbuf+AC.cbufnum;
1232  CBUF *CC = cbuf+AT.ebufnum;
1233  WORD *oldworkpointer = AT.WorkPointer, *termout;
1234 
1235  termout = AT.WorkPointer;
1236  NewSort(BHEAD0);
1237  while ( *in ) {
1238  if ( ConvertFromPoly(BHEAD in,termout,numxsymbol,CC->numrhs-startebuf+numxsymbol,startebuf-numxsymbol,1) <= 0 ) {
1239  LowerSortLevel();
1240  goto CalledFrom;
1241  }
1242  in += *in;
1243  AT.WorkPointer = termout + *termout;
1244 /*
1245  ConvertFromPoly leaves terms with subexpressions. Hence:
1246 */
1247  if ( Generator(BHEAD termout,C->numlhs) ) {
1248  LowerSortLevel();
1249  goto CalledFrom;
1250  }
1251  }
1252  AT.WorkPointer = oldworkpointer;
1253  if ( EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 ) goto CalledFrom;
1254  return(termout);
1255 
1256 CalledFrom:
1257  MLOCK(ErrorMessageLock);
1258  MesCall("TakeExtraSymbols");
1259  MUNLOCK(ErrorMessageLock);
1260  return(0);
1261 }
1262 
1263 /*
1264  #] TakeExtraSymbols :
1265  #[ MultiplyWithTerm :
1266 */
1267 
1268 WORD *MultiplyWithTerm(PHEAD WORD *in, WORD *term)
1269 {
1270  WORD *termout, *t, *tt, *tstop, *ttstop;
1271  WORD length, length1, length2;
1272  termout = AT.WorkPointer;
1273  NewSort(BHEAD0);
1274  while ( *in ) {
1275  tt = termout + 1;
1276  tstop = in + *in; tstop -= ABS(tstop[-1]); t = in + 1;
1277  while ( t < tstop ) *tt++ = *t++;
1278  ttstop = term + *term; ttstop -= ABS(ttstop[-1]); t = term + 1;
1279  while ( t < ttstop ) *tt++ = *t++;
1280  length1 = REDLENG(in[*in-1]); length2 = REDLENG(term[*term-1]);
1281  if ( MulRat(BHEAD (UWORD *)tstop,length1,
1282  (UWORD *)ttstop,length2,(UWORD *)tt,&length) ) goto CalledFrom;
1283  length = INCLENG(length);
1284  tt += ABS(length); tt[-1] = length;
1285  *termout = tt - termout;
1286  Normalize(BHEAD termout);
1287  StoreTerm(BHEAD termout);
1288  in += *in;
1289  }
1290  if ( EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 ) goto CalledFrom;
1291  return(termout);
1292 
1293 CalledFrom:
1294  MLOCK(ErrorMessageLock);
1295  MesCall("MultiplyWithTerm");
1296  MUNLOCK(ErrorMessageLock);
1297  return(0);
1298 }
1299 
1300 /*
1301  #] MultiplyWithTerm :
1302  #[ TakeContent :
1303 */
1316 WORD *TakeContent(PHEAD WORD *in, WORD *term)
1317 {
1318  GETBIDENTITY
1319  WORD *t, *tstop, *tcom, *tout, *tstore, *r, *rstop, *m, *mm, *w, *ww, *wterm;
1320  WORD *tnext, *tt, *tterm, code[2];
1321  WORD *inp, a, *den;
1322  int i, j, k, action = 0, sign;
1323  UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMbuffer2, *ap;
1324  WORD GCDlen, GCDlen2, LCMlen, LCMlen2, length, redlength, len1, len2;
1325  tout = tstore = term+1;
1326 /*
1327  #[ INDEX :
1328 */
1329  t = in;
1330  tnext = t + *t;
1331  tstop = tnext-ABS(tnext[-1]);
1332  t++;
1333  while ( t < tstop ) {
1334  if ( *t == INDEX ) {
1335  i = t[1]; NCOPY(tout,t,i); break;
1336  }
1337  else t += t[1];
1338  }
1339  if ( tout > tstore ) { /* There are indices in the first term */
1340  t = tnext;
1341  while ( *t ) {
1342  tnext = t + *t;
1343  rstop = tnext - ABS(tnext[-1]);
1344  r = t+1;
1345  while ( r < rstop ) {
1346  if ( *r != INDEX ) { r += r[1]; continue; }
1347  m = tstore+2;
1348  while ( m < tout ) {
1349  for ( i = 2; i < r[1]; i++ ) {
1350  if ( *m == r[i] ) break;
1351  if ( *m > r[i] ) continue;
1352  mm = m+1;
1353  while ( mm < tout ) { mm[-1] = mm[0]; mm++; }
1354  tout--; tstore[1]--; m--;
1355  break;
1356  }
1357  m++;
1358  }
1359  }
1360  if ( r >= rstop || tout <= tstore+2 ) {
1361  tout = tstore; break;
1362  }
1363  }
1364  if ( tout > tstore+2 ) { /* Now we have to take out what is in tstore */
1365  t = in; w = in;
1366  while ( *t ) {
1367  wterm = w;
1368  tnext = t + *t; t++; w++;
1369  while ( *t != INDEX ) { i = t[1]; NCOPY(w,t,i); }
1370  tt = t + t[1]; t += 2; r = tstore+2; ww = w; *w++ = INDEX; w++;
1371  while ( r < tout && t < tt ) {
1372  if ( *r > *t ) { *w++ = *t++; }
1373  else if ( *r == *t ) { r++; t++; }
1374  else goto CalledFrom;
1375  }
1376  if ( r < tout ) goto CalledFrom;
1377  while ( t < tt ) *w++ = *t++;
1378  ww[1] = w - ww;
1379  if ( ww[1] == 2 ) w = ww;
1380  while ( t < tnext ) *w++ = *t++;
1381  *wterm = w - wterm;
1382  }
1383  *w = 0;
1384  }
1385  tstore = tout;
1386  }
1387 /*
1388  #] INDEX :
1389  #[ VECTOR/DELTA :
1390 */
1391  code[0] = VECTOR; code[1] = DELTA;
1392  for ( k = 0; k < 2; k++ ) {
1393  t = in;
1394  tnext = t + *t;
1395  tstop = tnext-ABS(tnext[-1]);
1396  t++;
1397  while ( t < tstop ) {
1398  if ( *t == code[k] ) {
1399  i = t[1]; NCOPY(tout,t,i); break;
1400  }
1401  else t += t[1];
1402  }
1403  if ( tout > tstore ) { /* There are vectors in the first term */
1404  t = tnext;
1405  while ( *t ) {
1406  tnext = t + *t;
1407  rstop = tnext - ABS(tnext[-1]);
1408  r = t+1;
1409  while ( r < rstop ) {
1410  if ( *r != code[k] ) { r += r[1]; continue; }
1411  m = tstore+2;
1412  while ( m < tout ) {
1413  for ( i = 2; i < r[1]; i += 2 ) {
1414  if ( *m == r[i] && m[1] == r[i+1] ) break;
1415  if ( *m > r[i] || ( *m == r[i] && m[1] > r[i+1] ) ) continue;
1416  mm = m+2;
1417  while ( mm < tout ) { mm[-2] = mm[0]; mm[-1] = mm[1]; mm += 2; }
1418  tout -= 2; tstore[1] -= 2; m -= 2;
1419  break;
1420  }
1421  m += 2;
1422  }
1423  }
1424  if ( r >= rstop || tout <= tstore+2 ) {
1425  tout = tstore; break;
1426  }
1427  }
1428  if ( tout > tstore+2 ) { /* Now we have to take out what is in tstore */
1429  t = in; w = in;
1430  while ( *t ) {
1431  wterm = w;
1432  tnext = t + *t; t++; w++;
1433  while ( *t != code[k] ) { i = t[1]; NCOPY(w,t,i); }
1434  tt = t + t[1]; t += 2; r = tstore+2; ww = w; *w++ = code[k]; w++;
1435  while ( r < tout && t < tt ) {
1436  if ( ( *r > *t ) || ( *r == *t && r[1] > t[1] ) )
1437  { *w++ = *t++; *w++ = *t++; }
1438  else if ( *r == *t && r[1] == t[1] ) { r += 2; t += 2; }
1439  else goto CalledFrom;
1440  }
1441  if ( r < tout ) goto CalledFrom;
1442  while ( t < tt ) *w++ = *t++;
1443  ww[1] = w - ww;
1444  if ( ww[1] == 2 ) w = ww;
1445  while ( t < tnext ) *w++ = *t++;
1446  *wterm = w - wterm;
1447  }
1448  *w = 0;
1449  }
1450  tstore = tout;
1451  }
1452  }
1453 /*
1454  #] VECTOR/DELTA :
1455  #[ FUNCTIONS :
1456 */
1457  t = in;
1458  tnext = t + *t;
1459  tstop = tnext-ABS(tnext[-1]);
1460  t++;
1461  tcom = 0;
1462  while ( t < tstop ) {
1463  if ( *t >= FUNCTION ) {
1464  if ( functions[*t-FUNCTION].commute ) {
1465  if ( tcom == 0 ) { tcom = tstore; }
1466  else {
1467  for ( i = 0; i < t[1]; i++ ) {
1468  if ( t[i] != tcom[i] ) {
1469  MLOCK(ErrorMessageLock);
1470  MesPrint("GCD or factorization of more than one noncommuting object not allowed");
1471  MUNLOCK(ErrorMessageLock);
1472  goto CalledFrom;
1473  }
1474  }
1475  }
1476  }
1477  i = t[1]; NCOPY(tstore,t,i);
1478  }
1479  else t += t[1];
1480  }
1481  if ( tout > tstore ) {
1482  t = tnext;
1483  while ( *t ) {
1484  tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1485  r = tstore;
1486  while ( r < tout ) {
1487  tt = t;
1488  while ( tt < tstop ) {
1489  for ( i = 0; i < r[1]; i++ ) {
1490  if ( r[i] != tt[i] ) break;
1491  }
1492  if ( i == r[1] ) { r += r[1]; goto nextr1; }
1493  }
1494 /*
1495  Not encountered in this term. Scratch from list
1496 */
1497  m = r; mm = r + r[1];
1498  while ( mm < tout ) *m++ = *mm++;
1499  tout = m;
1500 nextr1:;
1501  }
1502  if ( tout <= tstore ) break;
1503  t += *t;
1504  }
1505  }
1506  if ( tout > tstore ) {
1507 /*
1508  Now we have one or more functions left that are common in all terms.
1509  Take them out. We do this one by one.
1510 */
1511  r = tstore;
1512  while ( r < tout ) {
1513  t = in; ww = in; w = ww+1;
1514  while ( *t ) {
1515  tnext = t + *t;
1516  t++;
1517  for(;;) {
1518  for ( i = 0; i < r[1]; i++ ) {
1519  if ( t[i] != r[i] ) {
1520  j = t[1]; NCOPY(w,t,j);
1521  break;
1522  }
1523  }
1524  if ( i == r[1] ) {
1525  t += t[1];
1526  while ( t < tnext ) *w++ = *t++;
1527  *ww = w - ww;
1528  break;
1529  }
1530  }
1531  }
1532  r += r[1];
1533  *w = 0;
1534  }
1535  tstore = tout;
1536  }
1537 /*
1538  #] FUNCTIONS :
1539  #[ SYMBOL :
1540 
1541  We make a list of symbols and their minimal powers.
1542  This includes negative powers. In the end we have to multiply by the
1543  inverse of this list. That takes out all negative powers and leaves
1544  things ready for further processing.
1545 */
1546  tterm = AT.WorkPointer; tt = tterm+1;
1547  tout[0] = SYMBOL; tout[1] = 2;
1548  t = in;
1549  while ( *t ) {
1550  tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1551  while ( t < tstop ) {
1552  if ( *t == SYMBOL ) {
1553  MergeSymbolLists(BHEAD tout,t,-1);
1554  break;
1555  }
1556  t += t[1];
1557  }
1558  t = tnext;
1559  }
1560  if ( tout[1] > 2 ) {
1561  t = tout;
1562  tt[0] = t[0]; tt[1] = t[1];
1563  for ( i = 2; i < t[1]; i += 2 ) {
1564  tt[i] = t[i]; tt[i+1] = -t[i+1];
1565  }
1566  tt += tt[1];
1567  tout += tout[1];
1568  action++;
1569  }
1570 /*
1571  #] SYMBOL :
1572  #[ DOTPRODUCT :
1573 
1574  We make a list of dotproducts and their minimal powers.
1575  This includes negative powers. In the end we have to multiply by the
1576  inverse of this list. That takes out all negative powers and leaves
1577  things ready for further processing.
1578 */
1579  tout[0] = DOTPRODUCT; tout[1] = 2;
1580  t = in;
1581  while ( *t ) {
1582  tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1583  while ( t < tstop ) {
1584  if ( *t == DOTPRODUCT ) {
1585  MergeDotproductLists(BHEAD tout,t,-1);
1586  break;
1587  }
1588  t += t[1];
1589  }
1590  t = tnext;
1591  }
1592  if ( tout[1] > 2 ) {
1593  t = tout;
1594  tt[0] = t[0]; tt[1] = t[1];
1595  for ( i = 2; i < t[1]; i += 3 ) {
1596  tt[i] = t[i]; tt[i+1] = t[i+1]; tt[i+2] = -t[i+2];
1597  }
1598  tt += tt[1];
1599  tout += tout[1];
1600  action++;
1601  }
1602 /*
1603  #] DOTPRODUCT :
1604  #[ Coefficient :
1605 
1606  Now we have to collect the GCD of the numerators
1607  and the LCM of the denominators.
1608 */
1609  AT.WorkPointer = tt;
1610  if ( AN.cmod != 0 ) {
1611  WORD x, ix, ip;
1612  t = in; tnext = t + *t; tstop = tnext - ABS(tnext[-1]);
1613  x = tstop[0];
1614  if ( tnext[-1] < 0 ) x += AC.cmod[0];
1615  if ( GetModInverses(x,(WORD)(AN.cmod[0]),&ix,&ip) ) goto CalledFrom;
1616  *tout++ = x; *tout++ = 1; *tout++ = 3;
1617  *tt++ = ix; *tt++ = 1; *tt++ = 3;
1618  }
1619  else {
1620  GCDbuffer = NumberMalloc("MakeInteger");
1621  GCDbuffer2 = NumberMalloc("MakeInteger");
1622  LCMbuffer = NumberMalloc("MakeInteger");
1623  LCMbuffer2 = NumberMalloc("MakeInteger");
1624  t = in;
1625  tnext = t + *t; length = tnext[-1];
1626  if ( length < 0 ) { sign = -1; length = -length; }
1627  else { sign = 1; }
1628  tstop = tnext - length;
1629  redlength = (length-1)/2;
1630  for ( i = 0; i < redlength; i++ ) {
1631  GCDbuffer[i] = (UWORD)(tstop[i]);
1632  LCMbuffer[i] = (UWORD)(tstop[redlength+i]);
1633  }
1634  GCDlen = LCMlen = redlength;
1635  while ( GCDbuffer[GCDlen-1] == 0 ) GCDlen--;
1636  while ( LCMbuffer[LCMlen-1] == 0 ) LCMlen--;
1637  t = tnext;
1638  while ( *t ) {
1639  tnext = t + *t; length = ABS(tnext[-1]);
1640  tstop = tnext - length; redlength = (length-1)/2;
1641  len1 = len2 = redlength;
1642  den = tstop + redlength;
1643  while ( tstop[len1-1] == 0 ) len1--;
1644  while ( den[len2-1] == 0 ) len2--;
1645  if ( GCDlen == 1 && GCDbuffer[0] == 1 ) {}
1646  else {
1647  GcdLong(BHEAD (UWORD *)tstop,len1,GCDbuffer,GCDlen,GCDbuffer2,&GCDlen2);
1648  ap = GCDbuffer; GCDbuffer = GCDbuffer2; GCDbuffer2 = ap;
1649  a = GCDlen; GCDlen = GCDlen2; GCDlen2 = a;
1650  }
1651  if ( len2 == 1 && den[0] == 1 ) {}
1652  else {
1653  GcdLong(BHEAD LCMbuffer,LCMlen,(UWORD *)den,len2,LCMbuffer2,&LCMlen2);
1654  DivLong((UWORD *)den,len2,LCMbuffer2,LCMlen2,
1655  GCDbuffer2,&GCDlen2,(UWORD *)AT.WorkPointer,&a);
1656  MulLong(LCMbuffer,LCMlen,GCDbuffer2,GCDlen2,LCMbuffer2,&LCMlen2);
1657  ap = LCMbuffer; LCMbuffer = LCMbuffer2; LCMbuffer2 = ap;
1658  a = LCMlen; LCMlen = LCMlen2; LCMlen2 = a;
1659  }
1660  t = tnext;
1661  }
1662  if ( GCDlen != 1 || GCDbuffer[0] != 1 || LCMlen != 1 || LCMbuffer[0] != 1 ) {
1663  redlength = GCDlen; if ( LCMlen > GCDlen ) redlength = LCMlen;
1664  for ( i = 0; i < GCDlen; i++ ) *tout++ = (WORD)(GCDbuffer[i]);
1665  for ( ; i < redlength; i++ ) *tout++ = 0;
1666  for ( i = 0; i < LCMlen; i++ ) *tout++ = (WORD)(LCMbuffer[i]);
1667  for ( ; i < redlength; i++ ) *tout++ = 0;
1668  *tout++ = (2*redlength+1)*sign;
1669  for ( i = 0; i < LCMlen; i++ ) *tt++ = (WORD)(LCMbuffer[i]);
1670  for ( ; i < redlength; i++ ) *tt++ = 0;
1671  for ( i = 0; i < GCDlen; i++ ) *tt++ = (WORD)(GCDbuffer[i]);
1672  for ( ; i < redlength; i++ ) *tt++ = 0;
1673  *tt++ = (2*redlength+1)*sign;
1674  action++;
1675  }
1676  else {
1677  *tout++ = 1; *tout++ = 1; *tout++ = 3*sign;
1678  *tt++ = 1; *tt++ = 1; *tt++ = 3*sign;
1679  if ( sign != 1 ) action++;
1680  }
1681  NumberFree(LCMbuffer2,"MakeInteger");
1682  NumberFree(LCMbuffer ,"MakeInteger");
1683  NumberFree(GCDbuffer2,"MakeInteger");
1684  NumberFree(GCDbuffer ,"MakeInteger");
1685  }
1686 /*
1687  #] Coefficient :
1688  #[ Multiply by the inverse content :
1689 */
1690  if ( action ) {
1691  *tterm = tt - tterm;
1692  AT.WorkPointer = tt;
1693  inp = MultiplyWithTerm(BHEAD in,tterm);
1694  AT.WorkPointer = tterm;
1695  in = inp;
1696  }
1697 /*
1698  #] Multiply by the inverse content :
1699 */
1700  *term = tout - term;
1701  AT.WorkPointer = tterm;
1702  return(in);
1703 CalledFrom:
1704  MLOCK(ErrorMessageLock);
1705  MesCall("TakeContent");
1706  MUNLOCK(ErrorMessageLock);
1707  return(0);
1708 }
1709 
1710 /*
1711  #] TakeContent :
1712  #[ MergeSymbolLists :
1713 
1714  Merges the extra list into the old.
1715  If par == -1 we take minimum powers
1716  If par == 1 we take maximum powers
1717  If par == 0 we take minimum of the absolute value of the powers
1718  if one is positive and the other negative we get zero.
1719  We assume that the symbols are in order in both lists
1720 */
1721 
1722 int MergeSymbolLists(PHEAD WORD *old, WORD *extra, int par)
1723 {
1724  GETBIDENTITY
1725  WORD *new = TermMalloc("MergeSymbolLists");
1726  WORD *t1, *t2, *fill;
1727  int i1,i2;
1728  fill = new + 2;
1729  i1 = old[1] - 2; i2 = extra[1] - 2;
1730  t1 = old + 2; t2 = extra + 2;
1731  switch ( par ) {
1732  case -1:
1733  while ( i1 > 0 && i2 > 0 ) {
1734  if ( *t1 > *t2 ) {
1735  if ( t2[1] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1736  else t2 += 2;
1737  }
1738  else if ( *t1 < *t2 ) {
1739  if ( t1[1] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1740  else t1 += 2;
1741  }
1742  else if ( t1[1] < t2[1] ) {
1743  *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1744  }
1745  else {
1746  *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1747  }
1748  }
1749  break;
1750  case 1:
1751  while ( i1 > 0 && i2 > 0 ) {
1752  if ( *t1 > *t2 ) {
1753  if ( t2[1] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1754  else t2 += 2;
1755  }
1756  else if ( *t1 < *t2 ) {
1757  if ( t1[1] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1758  else t1 += 2;
1759  }
1760  else if ( t1[1] > t2[1] ) {
1761  *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1762  }
1763  else {
1764  *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1765  }
1766  }
1767  break;
1768  case 0:
1769  while ( i1 > 0 && i2 > 0 ) {
1770  if ( *t1 > *t2 ) {
1771  t2 += 2;
1772  }
1773  else if ( *t1 < *t2 ) {
1774  t1 += 2;
1775  }
1776  else if ( ( t1[1] > 0 ) && ( t2[1] < 0 ) ) { t1 += 2; t2 += 2; }
1777  else if ( ( t1[1] < 0 ) && ( t2[1] > 0 ) ) { t1 += 2; t2 += 2; }
1778  else if ( t1[1] > 0 ) {
1779  if ( t1[1] < t2[1] ) {
1780  *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1781  }
1782  else {
1783  *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1784  }
1785  }
1786  else {
1787  if ( t2[1] < t1[1] ) {
1788  *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1789  }
1790  else {
1791  *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1792  }
1793  }
1794  }
1795  break;
1796  }
1797  i1 = new[1] = fill - new;
1798  t2 = new; t1 = old; NCOPY(t1,t2,i1);
1799  TermFree(new,"MergeSymbolLists");
1800  return(0);
1801 }
1802 
1803 /*
1804  #] MergeSymbolLists :
1805  #[ MergeDotproductLists :
1806 
1807  Merges the extra list into the old.
1808  If par == -1 we take minimum powers
1809  If par == 1 we take maximum powers
1810  If par == 0 we take minimum of the absolute value of the powers
1811  if one is positive and the other negative we get zero.
1812  We assume that the dotproducts are in order in both lists
1813 */
1814 
1815 int MergeDotproductLists(PHEAD WORD *old, WORD *extra, int par)
1816 {
1817  GETBIDENTITY
1818  WORD *new = TermMalloc("MergeDotproductLists");
1819  WORD *t1, *t2, *fill;
1820  int i1,i2;
1821  fill = new + 2;
1822  i1 = old[1] - 2; i2 = extra[1] - 2;
1823  t1 = old + 2; t2 = extra + 2;
1824  switch ( par ) {
1825  case -1:
1826  while ( i1 > 0 && i2 > 0 ) {
1827  if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1828  if ( t2[2] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; }
1829  else t2 += 3;
1830  }
1831  else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1832  if ( t1[2] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; }
1833  else t1 += 3;
1834  }
1835  else if ( t1[2] < t2[2] ) {
1836  *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1837  }
1838  else {
1839  *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1840  }
1841  }
1842  break;
1843  case 1:
1844  while ( i1 > 0 && i2 > 0 ) {
1845  if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1846  if ( t2[2] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; }
1847  else t2 += 3;
1848  }
1849  else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1850  if ( t1[2] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; }
1851  else t1 += 3;
1852  }
1853  else if ( t1[2] > t2[2] ) {
1854  *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1855  }
1856  else {
1857  *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1858  }
1859  }
1860  break;
1861  case 0:
1862  while ( i1 > 0 && i2 > 0 ) {
1863  if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1864  t2 += 3;
1865  }
1866  else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1867  t1 += 3;
1868  }
1869  else if ( ( t1[2] > 0 ) && ( t2[2] < 0 ) ) { t1 += 3; t2 += 3; }
1870  else if ( ( t1[2] < 0 ) && ( t2[2] > 0 ) ) { t1 += 3; t2 += 3; }
1871  else if ( t1[2] > 0 ) {
1872  if ( t1[2] < t2[2] ) {
1873  *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1874  }
1875  else {
1876  *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1877  }
1878  }
1879  else {
1880  if ( t2[2] < t1[2] ) {
1881  *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1882  }
1883  else {
1884  *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1885  }
1886  }
1887  }
1888  break;
1889  }
1890  i1 = new[1] = fill - new;
1891  t2 = new; t1 = old; NCOPY(t1,t2,i1);
1892  TermFree(new,"MergeDotproductLists");
1893  return(0);
1894 }
1895 
1896 /*
1897  #] MergeDotproductLists :
1898  #[ CreateExpression :
1899 
1900  Looks for the expression in the argument, reads it and puts it
1901  in a buffer. Returns the address of the buffer.
1902  We send the expression through the Generator system, because there
1903  may be unsubstituted (sub)expressions as in
1904  Local F = (a+b);
1905  Local G = gcd_(F,...);
1906 */
1907 
1908 WORD *CreateExpression(PHEAD WORD nexp)
1909 {
1910  GETBIDENTITY
1911  CBUF *C = cbuf+AC.cbufnum;
1912  POSITION startposition, oldposition;
1913  FILEHANDLE *fi;
1914  WORD *term, *oldipointer = AR.CompressPointer;
1915 ;
1916  switch ( Expressions[nexp].status ) {
1917  case HIDDENLEXPRESSION:
1918  case HIDDENGEXPRESSION:
1919  case DROPHLEXPRESSION:
1920  case DROPHGEXPRESSION:
1921  case UNHIDELEXPRESSION:
1922  case UNHIDEGEXPRESSION:
1923  AR.GetOneFile = 2; fi = AR.hidefile;
1924  break;
1925  default:
1926  AR.GetOneFile = 0; fi = AR.infile;
1927  break;
1928  }
1929  SeekScratch(fi,&oldposition);
1930  startposition = AS.OldOnFile[nexp];
1931  term = AT.WorkPointer;
1932  if ( GetOneTerm(BHEAD term,fi,&startposition,0) <= 0 ) goto CalledFrom;
1933  NewSort(BHEAD0);
1934  AR.CompressPointer = oldipointer;
1935  while ( GetOneTerm(BHEAD term,fi,&startposition,0) > 0 ) {
1936  AT.WorkPointer = term + *term;
1937  if ( Generator(BHEAD term,C->numlhs) ) {
1938  LowerSortLevel();
1939  goto CalledFrom;
1940  }
1941  AR.CompressPointer = oldipointer;
1942  }
1943  AT.WorkPointer = term;
1944  if ( EndSort(BHEAD (WORD *)((VOID *)(&term)),2) < 0 ) goto CalledFrom;
1945  SetScratch(fi,&oldposition);
1946  return(term);
1947 CalledFrom:
1948  MLOCK(ErrorMessageLock);
1949  MesCall("CreateExpression");
1950  MUNLOCK(ErrorMessageLock);
1951  Terminate(-1);
1952  return(0);
1953 }
1954 
1955 /*
1956  #] CreateExpression :
1957  #[ GCDterms : GCD of two terms
1958 
1959  Computes the GCD of two terms.
1960  Output in termout.
1961  termout may overlap with term1.
1962 */
1963 
1964 int GCDterms(PHEAD WORD *term1, WORD *term2, WORD *termout)
1965 {
1966  GETBIDENTITY
1967  WORD *t1, *t1stop, *t1next, *t2, *t2stop, *t2next, *tout, *tt1, *tt2;
1968  int count1, count2, i, ii, x1, sign;
1969  WORD length1, length2;
1970  t1 = term1 + *term1; t1stop = t1 - ABS(t1[-1]); t1 = term1+1;
1971  t2 = term2 + *term2; t2stop = t2 - ABS(t2[-1]); t2 = term2+1;
1972  tout = termout+1;
1973  while ( t1 < t1stop ) {
1974  t1next = t1 + t1[1];
1975  t2 = term2+1;
1976  if ( *t1 == SYMBOL ) {
1977  while ( t2 < t2stop && *t2 != SYMBOL ) t2 += t2[1];
1978  if ( *t2 == SYMBOL ) {
1979  t2next = t2+t2[1];
1980  tt1 = t1+2; tt2 = t2+2; count1 = 0;
1981  while ( tt1 < t1next && tt2 < t2next ) {
1982  if ( *tt1 < *tt2 ) tt1 += 2;
1983  else if ( *tt1 > *tt2 ) tt2 += 2;
1984  else if ( ( tt1[1] > 0 && tt2[1] < 0 ) ||
1985  ( tt2[1] > 0 && tt1[1] < 0 ) ) {
1986  tt1 += 2; tt2 += 2;
1987  }
1988  else {
1989  x1 = tt1[1];
1990  if ( tt1[1] < 0 ) { if ( tt2[1] > x1 ) x1 = tt2[1]; }
1991  else { if ( tt2[1] < x1 ) x1 = tt2[1]; }
1992  tout[count1+2] = *tt1;
1993  tout[count1+3] = x1;
1994  tt1 += 2; tt2 += 2;
1995  count1 += 2;
1996  }
1997  }
1998  if ( count1 > 0 ) {
1999  *tout = SYMBOL; tout[1] = count1+2; tout += tout[1];
2000  }
2001  }
2002  }
2003  else if ( *t1 == DOTPRODUCT ) {
2004  while ( t2 < t2stop && *t2 != DOTPRODUCT ) t2 += t2[1];
2005  if ( *t2 == DOTPRODUCT ) {
2006  t2next = t2+t2[1];
2007  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2008  while ( tt1 < t1next && tt2 < t2next ) {
2009  if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 3;
2010  else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 3;
2011  else if ( ( tt1[2] > 0 && tt2[2] < 0 ) ||
2012  ( tt2[2] > 0 && tt1[2] < 0 ) ) {
2013  tt1 += 3; tt2 += 3;
2014  }
2015  else {
2016  x1 = tt1[2];
2017  if ( tt1[2] < 0 ) { if ( tt2[2] > x1 ) x1 = tt2[2]; }
2018  else { if ( tt2[2] < x1 ) x1 = tt2[2]; }
2019  tout[count1+2] = *tt1;
2020  tout[count1+3] = tt1[1];
2021  tout[count1+4] = x1;
2022  tt1 += 3; tt2 += 3;
2023  count1 += 3;
2024  }
2025  }
2026  if ( count1 > 0 ) {
2027  *tout = DOTPRODUCT; tout[1] = count1+2; tout += tout[1];
2028  }
2029  }
2030  }
2031  else if ( *t1 == VECTOR ) {
2032  while ( t2 < t2stop && *t2 != VECTOR ) t2 += t2[1];
2033  if ( *t2 == VECTOR ) {
2034  t2next = t2+t2[1];
2035  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2036  while ( tt1 < t1next && tt2 < t2next ) {
2037  if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2038  else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2039  else {
2040  tout[count1+2] = *tt1;
2041  tout[count1+3] = tt1[1];
2042  tt1 += 2; tt2 += 2;
2043  count1 += 2;
2044  }
2045  }
2046  if ( count1 > 0 ) {
2047  *tout = VECTOR; tout[1] = count1+2; tout += tout[1];
2048  }
2049  }
2050  }
2051  else if ( *t1 == INDEX ) {
2052  while ( t2 < t2stop && *t2 != INDEX ) t2 += t2[1];
2053  if ( *t2 == INDEX ) {
2054  t2next = t2+t2[1];
2055  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2056  while ( tt1 < t1next && tt2 < t2next ) {
2057  if ( *tt1 < *tt2 ) tt1 += 1;
2058  else if ( *tt1 > *tt2 ) tt2 += 1;
2059  else {
2060  tout[count1+2] = *tt1;
2061  tt1 += 1; tt2 += 1;
2062  count1 += 1;
2063  }
2064  }
2065  if ( count1 > 0 ) {
2066  *tout = INDEX; tout[1] = count1+2; tout += tout[1];
2067  }
2068  }
2069  }
2070  else if ( *t1 == DELTA ) {
2071  while ( t2 < t2stop && *t2 != DELTA ) t2 += t2[1];
2072  if ( *t2 == DELTA ) {
2073  t2next = t2+t2[1];
2074  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2075  while ( tt1 < t1next && tt2 < t2next ) {
2076  if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2077  else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2078  else {
2079  tout[count1+2] = *tt1;
2080  tout[count1+3] = tt1[1];
2081  tt1 += 2; tt2 += 2;
2082  count1 += 2;
2083  }
2084  }
2085  if ( count1 > 0 ) {
2086  *tout = DELTA; tout[1] = count1+2; tout += tout[1];
2087  }
2088  }
2089  }
2090  else if ( *t1 >= FUNCTION ) { /* noncommuting functions? Forbidden! */
2091 /*
2092  Count how many times this function occurs.
2093  Then count how many times it is in term2.
2094 */
2095  count1 = 1;
2096  while ( t1next < t1stop && *t1 == *t1next && t1[1] == t1next[1] ) {
2097  for ( i = 2; i < t1[1]; i++ ) {
2098  if ( t1[i] != t1next[i] ) break;
2099  }
2100  if ( i < t1[1] ) break;
2101  count1++;
2102  t1next += t1next[1];
2103  }
2104  count2 = 0;
2105  while ( t2 < t2stop ) {
2106  if ( *t2 == *t1 && t2[1] == t1[1] ) {
2107  for ( i = 2; i < t1[1]; i++ ) {
2108  if ( t2[i] != t1[i] ) break;
2109  }
2110  if ( i >= t1[1] ) count2++;
2111  }
2112  t2 += t2[1];
2113  }
2114  if ( count1 < count2 ) count2 = count1; /* number of common occurrences */
2115  if ( count2 > 0 ) {
2116  if ( tout == t1 ) {
2117  while ( count2 > 0 ) { tout += tout[1]; count2--; }
2118  }
2119  else {
2120  i = t1[1]*count2;
2121  NCOPY(tout,t1,i);
2122  }
2123  }
2124  }
2125  t1 = t1next;
2126  }
2127 /*
2128  Now the coefficients. They are in t1stop and t2stop. Should go to tout.
2129 */
2130  sign = 1;
2131  length1 = term1[*term1-1]; ii = i = ABS(length1); t1 = t1stop;
2132  if ( t1 != tout ) { NCOPY(tout,t1,i); tout -= ii; }
2133  length2 = term2[*term2-1];
2134  if ( length1 < 0 && length2 < 0 ) sign = -1;
2135  if ( AccumGCD(BHEAD (UWORD *)tout,&length1,(UWORD *)t2stop,length2) ) {
2136  MLOCK(ErrorMessageLock);
2137  MesCall("GCDterms");
2138  MUNLOCK(ErrorMessageLock);
2139  SETERROR(-1)
2140  }
2141  if ( sign < 0 && length1 > 0 ) length1 = -length1;
2142  tout += ABS(length1); tout[-1] = length1;
2143  *termout = tout - termout;
2144  return(0);
2145 }
2146 
2147 /*
2148  #] GCDterms :
2149  #] GCDfunction :
2150  #[ DIVfunction :
2151 
2152  Input: a div_ function that has two arguments inside a term.
2153  Action: Calculates [arg1/arg2] using polynomial techniques if needed.
2154  Output: The output result is combined with the remainder of the term
2155  and sent to Generator for further processing.
2156  Note that the output can be just a number or many terms.
2157  In case par == 0 the output is [arg1/arg2]
2158  In case par == 1 the output is [arg1%arg2]
2159 */
2160 
2161 WORD divrem[3] = { DIVFUNCTION, REMFUNCTION, INVERSEFUNCTION };
2162 
2163 int DIVfunction(PHEAD WORD *term,WORD level,int par)
2164 {
2165  GETBIDENTITY
2166  WORD *t, *tt, *r, *arg1 = 0, *arg2 = 0, *arg3 = 0, *termout;
2167  WORD *tstop, *tend, *r3, *rr, *rstop, tlength, rlength, newlength;
2168  WORD *proper1, *proper2, *proper3 = 0;
2169  int numargs = 0, type1, type2, actionflag1, actionflag2;
2170  WORD startebuf = cbuf[AT.ebufnum].numrhs;
2171  if ( par < 0 || par > 2 ) {
2172  MLOCK(ErrorMessageLock);
2173  MesPrint("Internal error. Illegal parameter %d in DIVfunction.",par);
2174  MUNLOCK(ErrorMessageLock);
2175  Terminate(-1);
2176  }
2177 /*
2178  Find the function
2179 */
2180  tend = term + *term; tstop = tend - ABS(tend[-1]);
2181  t = term+1;
2182  while ( t < tstop ) {
2183  if ( *t != divrem[par] ) { t += t[1]; continue; }
2184  r = t + FUNHEAD;
2185  tt = t + t[1]; numargs = 0;
2186  while ( r < tt ) {
2187  if ( numargs == 0 ) { arg1 = r; }
2188  if ( numargs == 1 ) { arg2 = r; }
2189  numargs++;
2190  NEXTARG(r);
2191  }
2192  if ( numargs == 2 ) break;
2193  t = tt;
2194  }
2195  if ( t >= tstop ) {
2196  MLOCK(ErrorMessageLock);
2197  MesPrint("Internal error. Indicated div_ or rem_ function not encountered.");
2198  MUNLOCK(ErrorMessageLock);
2199  Terminate(-1);
2200  }
2201 /*
2202  We have two arguments in arg1 and arg2.
2203 */
2204  if ( *arg1 == -SNUMBER && arg1[1] == 0 ) {
2205  if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2206 zerozero:;
2207  MLOCK(ErrorMessageLock);
2208  MesPrint("0/0 in either div_ or rem_ function.");
2209  MUNLOCK(ErrorMessageLock);
2210  Terminate(-1);
2211  }
2212  return(0);
2213  }
2214  if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2215 divzero:;
2216  MLOCK(ErrorMessageLock);
2217  MesPrint("Division by zero in either div_ or rem_ function.");
2218  MUNLOCK(ErrorMessageLock);
2219  Terminate(-1);
2220  }
2221  if ( ( arg1 = ConvertArgument(BHEAD arg1, &type1) ) == 0 ) goto CalledFrom;
2222  if ( ( arg2 = ConvertArgument(BHEAD arg2, &type2) ) == 0 ) goto CalledFrom;
2223  if ( *arg1 == 0 ) {
2224  if ( *arg2 == 0 ) {
2225  M_free(arg2,"DIVfunction");
2226  M_free(arg1,"DIVfunction");
2227  goto zerozero;
2228  }
2229  M_free(arg2,"DIVfunction");
2230  M_free(arg1,"DIVfunction");
2231  return(0);
2232  }
2233  if ( *arg2 == 0 ) {
2234  M_free(arg2,"DIVfunction");
2235  M_free(arg1,"DIVfunction");
2236  goto divzero;
2237  }
2238  if ( ( proper1 = PutExtraSymbols(BHEAD arg1,startebuf,&actionflag1) ) == 0 ) goto CalledFrom;
2239  if ( ( proper2 = PutExtraSymbols(BHEAD arg2,startebuf,&actionflag2) ) == 0 ) goto CalledFrom;
2240 /*
2241  if ( type2 == 0 ) M_free(arg2,"DIVfunction");
2242  else {
2243  DOLLARS d = ((DOLLARS)arg2)-1;
2244  if ( d->factors ) M_free(d->factors,"Dollar factors");
2245  M_free(d,"Copy of dollar variable");
2246  }
2247 */
2248  M_free(arg2,"DIVfunction");
2249 /*
2250  if ( type1 == 0 ) M_free(arg1,"DIVfunction");
2251  else {
2252  DOLLARS d = ((DOLLARS)arg1)-1;
2253  if ( d->factors ) M_free(d->factors,"Dollar factors");
2254  M_free(d,"Copy of dollar variable");
2255  }
2256 */
2257  M_free(arg1,"DIVfunction");
2258  if ( par == 0 ) proper3 = poly_div(BHEAD proper1, proper2);
2259  else if ( par == 1 ) proper3 = poly_rem(BHEAD proper1, proper2);
2260  else if ( par == 2 ) proper3 = poly_inverse(BHEAD proper1, proper2);
2261  if ( proper3 == 0 ) goto CalledFrom;
2262  if ( actionflag1 || actionflag2 ) {
2263  if ( ( arg3 = TakeExtraSymbols(BHEAD proper3,startebuf) ) == 0 ) goto CalledFrom;
2264  M_free(proper3,"DIVfunction");
2265  }
2266  else {
2267  arg3 = proper3;
2268  }
2269  M_free(proper2,"DIVfunction");
2270  M_free(proper1,"DIVfunction");
2271  cbuf[AT.ebufnum].numrhs = startebuf;
2272  if ( *arg3 ) {
2273  termout = AT.WorkPointer;
2274  tlength = tend[-1];
2275  tlength = REDLENG(tlength);
2276  r3 = arg3;
2277  while ( *r3 ) {
2278  tt = term + 1; rr = termout + 1;
2279  while ( tt < t ) *rr++ = *tt++;
2280  r = r3 + 1;
2281  r3 = r3 + *r3;
2282  rstop = r3 - ABS(r3[-1]);
2283  while ( r < rstop ) *rr++ = *r++;
2284  tt += t[1];
2285  while ( tt < tstop ) *rr++ = *tt++;
2286  rlength = r3[-1];
2287  rlength = REDLENG(rlength);
2288  if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)rstop,rlength,
2289  (UWORD *)rr,&newlength) < 0 ) goto CalledFrom;
2290  rlength = INCLENG(newlength);
2291  rr += ABS(rlength);
2292  rr[-1] = rlength;
2293  *termout = rr - termout;
2294  AT.WorkPointer = rr;
2295  if ( Generator(BHEAD termout,level) ) goto CalledFrom;
2296  }
2297  AT.WorkPointer = termout;
2298  }
2299  M_free(arg3,"DIVfunction");
2300  return(0);
2301 CalledFrom:
2302  MLOCK(ErrorMessageLock);
2303  MesCall("DIVfunction");
2304  MUNLOCK(ErrorMessageLock);
2305  SETERROR(-1)
2306 }
2307 
2308 /*
2309  #] DIVfunction :
2310  #[ ConvertArgument :
2311 
2312  Converts an argument to a general notation in allocated space.
2313 */
2314 
2315 WORD *ConvertArgument(PHEAD WORD *arg, int *type)
2316 {
2317  WORD *output, *t, *r;
2318  int i, size;
2319  if ( *arg > 0 ) {
2320  output = (WORD *)Malloc1((*arg)*sizeof(WORD),"ConvertArgument");
2321  i = *arg - ARGHEAD; t = arg + ARGHEAD; r = output;
2322  NCOPY(r,t,i);
2323  *r = 0; *type = 0;
2324  return(output);
2325  }
2326  if ( *arg == -EXPRESSION ) {
2327  *type = 0;
2328  return(CreateExpression(BHEAD arg[1]));
2329  }
2330  if ( *arg == -DOLLAREXPRESSION ) {
2331  DOLLARS d;
2332  *type = 1;
2333  d = DolToTerms(BHEAD arg[1]);
2334 /*
2335  The problem is that DolToTerms creates a copy of the dollar variable.
2336  If we just return d->where we create a memory leak. Hence we have to
2337  copy the contents of d->where to a new buffer
2338 */
2339  output = (WORD *)Malloc1((d->size+1)*sizeof(WORD),"Copy of dollar content");
2340  WCOPY(output,d->where,d->size+1);
2341  if ( d->factors ) { M_free(d->factors,"Dollar factors"); d->factors = 0; }
2342  M_free(d,"Copy of dollar variable");
2343  return(output);
2344  }
2345 #if ( FUNHEAD > 4 )
2346  size = FUNHEAD+5;
2347 #else
2348  size = 9;
2349 #endif
2350  output = (WORD *)Malloc1(size*sizeof(WORD),"ConvertArgument");
2351  switch(*arg) {
2352  case -SYMBOL:
2353  output[0] = 8; output[1] = SYMBOL; output[2] = 4; output[3] = arg[1];
2354  output[4] = 1; output[5] = 1; output[6] = 1; output[7] = 3;
2355  output[8] = 0;
2356  break;
2357  case -INDEX:
2358  case -VECTOR:
2359  case -MINVECTOR:
2360  output[0] = 7; output[1] = INDEX; output[2] = 3; output[3] = arg[1];
2361  output[4] = 1; output[5] = 1;
2362  if ( *arg == -MINVECTOR ) output[6] = -3;
2363  else output[6] = 3;
2364  output[7] = 0;
2365  break;
2366  case -SNUMBER:
2367  output[0] = 4;
2368  if ( arg[1] < 0 ) {
2369  output[1] = -arg[1]; output[2] = 1; output[3] = -3;
2370  }
2371  else {
2372  output[1] = arg[1]; output[2] = 1; output[3] = 3;
2373  }
2374  output[4] = 0;
2375  break;
2376  default:
2377  output[0] = FUNHEAD+4;
2378  output[1] = -*arg;
2379  output[2] = FUNHEAD;
2380  for ( i = 3; i <= FUNHEAD; i++ ) output[i] = 0;
2381  output[FUNHEAD+1] = 1;
2382  output[FUNHEAD+2] = 1;
2383  output[FUNHEAD+3] = 3;
2384  output[FUNHEAD+4] = 0;
2385  break;
2386  }
2387  *type = 0;
2388  return(output);
2389 }
2390 
2391 /*
2392  #] ConvertArgument :
2393 
2394 */
2395 
2396 
2397 
WORD * TakeContent(PHEAD WORD *in, WORD *term)
Definition: ratio.c:1316
Definition: structs.h:618
#define PHEAD
Definition: ftypes.h:56
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
Definition: notation.c:508
Definition: structs.h:908
WORD InsertTerm(PHEAD WORD *, WORD, WORD, WORD *, WORD *, WORD)
Definition: proces.c:2350
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition: reken.c:1443
WORD StoreTerm(PHEAD WORD *)
Definition: sort.c:4070
WORD ** rhs
Definition: structs.h:913
VOID LowerSortLevel()
Definition: sort.c:4435
WORD * Buffer
Definition: structs.h:909
WORD NewSort(PHEAD0)
Definition: sort.c:553
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:2865
LONG EndSort(PHEAD WORD *, int)
Definition: sort.c:632