FORM  4.2.1
ratio.c
Go to the documentation of this file.
1 
8 /* #[ License : */
9 /*
10  * Copyright (C) 1984-2017 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))*2;
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, argsdone = 0, *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  argsdone++;
679  continue;
680  }
681  else if ( *tf != firstshort ) {
682  if ( *tf != -INDEX && *tf != -VECTOR && *tf != -MINVECTOR ) {
683  argsdone++; gcdisone = 1; break;
684  }
685  if ( firstshort != -INDEX && firstshort != -VECTOR && firstshort != -MINVECTOR ) {
686  argsdone++; gcdisone = 1; break;
687  }
688  if ( tf[1] != firstvalue ) {
689  argsdone++; gcdisone = 1; break;
690  }
691  if ( *t == -MINVECTOR ) { firstshort = -VECTOR; }
692  if ( firstshort == -MINVECTOR ) { firstshort = -VECTOR; }
693  }
694  else if ( *tf > -FUNCTION && *tf != -SNUMBER && tf[1] != firstvalue ) {
695  argsdone++; gcdisone = 1; 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  argsdone++;
711  if ( firstvalue == 1 ) { gcdisone = 1; break; }
712  }
713  NEXTARG(tf);
714  }
715  termout = AT.WorkPointer;
716  AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
717  if ( AT.WorkPointer > AT.WorkTop ) {
718  MLOCK(ErrorMessageLock);
719  MesWork();
720  MUNLOCK(ErrorMessageLock);
721  return(-1);
722  }
723 /*
724  #] Do short arguments :
725  #[ Do trivial GCD :
726 
727  Copy head
728 */
729  i = t - term; tin = term; tout = termout;
730  NCOPY(tout,tin,i);
731  if ( gcdisone || ( firstshort == -SNUMBER && firstvalue == 1 ) ) {
732  sign = 1;
733 gcdone:
734  tin += t[1]; tstop = term + *term;
735  while ( tin < tstop ) *tout++ = *tin++;
736  *termout = tout - termout;
737  if ( sign < 0 ) tout[-1] = -tout[-1];
738  AT.WorkPointer = tout;
739  if ( argsdone && Generator(BHEAD termout,level) < 0 ) goto CalledFrom;
740  AT.WorkPointer = termout;
741  AT.pWorkPointer = args;
742  return(0);
743  }
744 /*
745  #] Do trivial GCD :
746  #[ Do short argument GCD :
747 */
748  if ( numargs == 0 ) { /* basically we are done */
749 doshort:
750  sign = 1;
751  if ( firstshort == 0 ) goto gcdone;
752  if ( firstshort == -SNUMBER ) {
753  *tout++ = SNUMBER; *tout++ = 4; *tout++ = firstvalue; *tout++ = 1;
754  goto gcdone;
755  }
756  else if ( firstshort == -SYMBOL ) {
757  *tout++ = SYMBOL; *tout++ = 4; *tout++ = firstvalue; *tout++ = 1;
758  goto gcdone;
759  }
760  else if ( firstshort == -VECTOR || firstshort == -INDEX ) {
761  *tout++ = INDEX; *tout++ = 3; *tout++ = firstvalue; goto gcdone;
762  }
763  else if ( firstshort == -MINVECTOR ) {
764  sign = -1;
765  *tout++ = INDEX; *tout++ = 3; *tout++ = firstvalue; goto gcdone;
766  }
767  else if ( firstshort <= -FUNCTION ) {
768  *tout++ = firstvalue; *tout++ = FUNHEAD; FILLFUN(tout);
769  goto gcdone;
770  }
771  else {
772  MLOCK(ErrorMessageLock);
773  MesPrint("Internal error. Illegal short argument in GCDfunction.");
774  MUNLOCK(ErrorMessageLock);
775  Terminate(-1);
776  }
777  }
778 /*
779  #] Do short argument GCD :
780  #[ Convert short argument :
781 
782  Now we allocate space for the arguments in general notation.
783  First the special one if there were short arguments
784 */
785  if ( firstshort ) {
786  switch ( firstshort ) {
787  case -SNUMBER:
788  sh[0] = 4; sh[1] = ABS(firstvalue); sh[2] = 1;
789  if ( firstvalue < 0 ) sh[3] = -3;
790  else sh[3] = 3;
791  sh[4] = 0;
792  break;
793  case -MINVECTOR:
794  case -VECTOR:
795  case -INDEX:
796  sh[0] = 8; sh[1] = INDEX; sh[2] = 3; sh[3] = firstvalue;
797  sh[4] = 1; sh[5] = 1;
798  if ( firstshort == -MINVECTOR ) sh[6] = -3;
799  else sh[6] = 3;
800  sh[7] = 0;
801  break;
802  case -SYMBOL:
803  sh[0] = 8; sh[1] = SYMBOL; sh[2] = 4; sh[3] = firstvalue; sh[4] = 1;
804  sh[5] = 1; sh[6] = 1; sh[7] = 3; sh[8] = 0;
805  break;
806  default:
807  sh[0] = FUNHEAD+4; sh[1] = firstshort; sh[2] = FUNHEAD;
808  for ( i = 2; i < FUNHEAD; i++ ) sh[i+1] = 0;
809  sh[FUNHEAD+1] = 1; sh[FUNHEAD+2] = 1; sh[FUNHEAD+3] = 3; sh[FUNHEAD+4] = 0;
810  break;
811  }
812  }
813 /*
814  #] Convert short argument :
815  #[ Sort arguments :
816 
817  Now we should sort the arguments in a way that the dollars and the
818  expressions come last. That way we may never need them.
819 */
820  for ( i = 1; i < numargs; i++ ) {
821  for ( ii = i; ii > 0; ii-- ) {
822  arg1 = AT.pWorkSpace[args+ii];
823  arg2 = AT.pWorkSpace[args+ii-1];
824  if ( *arg1 < 0 ) {
825  if ( *arg2 < 0 ) {
826  if ( *arg1 == -EXPRESSION ) break;
827  if ( *arg2 == -DOLLAREXPRESSION ) break;
828  AT.pWorkSpace[args+ii] = arg2;
829  AT.pWorkSpace[args+ii-1] = arg1;
830  }
831  else break;
832  }
833  else if ( *arg2 < 0 ) {
834  AT.pWorkSpace[args+ii] = arg2;
835  AT.pWorkSpace[args+ii-1] = arg1;
836  }
837  else {
838  if ( *arg1 > *arg2 ) {
839  AT.pWorkSpace[args+ii] = arg2;
840  AT.pWorkSpace[args+ii-1] = arg1;
841  }
842  else break;
843  }
844  }
845  }
846 /*
847  #] Sort arguments :
848  #[ There is a single term argument :
849 */
850  if ( firstshort ) {
851  mh = sh; istart = 0;
852 oneterm:;
853  for ( i = istart; i < numargs; i++ ) {
854  arg1 = AT.pWorkSpace[args+i];
855  if ( *arg1 > 0 ) {
856  oldval1 = arg1[*arg1]; arg1[*arg1] = 0;
857  m = arg1+ARGHEAD;
858  while ( *m ) {
859  GCDterms(BHEAD mh,m,mh); m += *m;
860  if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
861  gcdisone = 1; sign = 1; arg1[*arg1] = oldval1; goto gcdone;
862  }
863  }
864  arg1[*arg1] = oldval1;
865  }
866  else if ( *arg1 == -DOLLAREXPRESSION ) {
867  if ( ( d = DolToTerms(BHEAD arg1[1]) ) != 0 ) {
868  m = d->where;
869  while ( *m ) {
870  GCDterms(BHEAD mh,m,mh); m += *m;
871  argsdone++;
872  if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
873  gcdisone = 1; sign = 1;
874  if ( d->factors ) M_free(d->factors,"Dollar factors");
875  M_free(d,"Copy of dollar variable"); goto gcdone;
876  }
877  }
878  if ( d->factors ) M_free(d->factors,"Dollar factors");
879  M_free(d,"Copy of dollar variable");
880  }
881  }
882  else {
883  mm = CreateExpression(BHEAD arg1[1]);
884  m = mm;
885  while ( *m ) {
886  GCDterms(BHEAD mh,m,mh); m += *m;
887  argsdone++;
888  if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
889  gcdisone = 1; sign = 1; M_free(mm,"CreateExpression"); goto gcdone;
890  }
891  }
892  M_free(mm,"CreateExpression");
893  }
894  }
895  if ( firstshort ) {
896  if ( mh[0] == 4 ) {
897  firstshort = -SNUMBER; firstvalue = mh[1] * (mh[3]/3);
898  }
899  else if ( mh[1] == SYMBOL ) {
900  firstshort = -SYMBOL; firstvalue = mh[3];
901  }
902  else if ( mh[1] == INDEX ) {
903  firstshort = -INDEX; firstvalue = mh[3];
904  if ( mh[6] == -3 ) firstshort = -MINVECTOR;
905  }
906  else if ( mh[1] >= FUNCTION ) {
907  firstshort = -mh[1]; firstvalue = mh[1];
908  }
909  goto doshort;
910  }
911  else {
912 /*
913  We have a GCD that is only a single term.
914  Paste it in and combine the coefficients.
915 */
916  mh[mh[0]] = 0;
917  mm = mh;
918  ii = 0;
919  goto multiterms;
920  }
921  }
922 /*
923  Now we have only regular arguments.
924  But some have not yet been expanded.
925  Check whether there are proper long arguments and if so if there is
926  one with just a single term
927 */
928  for ( i = 0; i < numargs; i++ ) {
929  arg1 = AT.pWorkSpace[args+i];
930  if ( *arg1 > 0 && arg1[ARGHEAD]+ARGHEAD == *arg1 ) {
931 /*
932  We have an argument with a single term
933 */
934  if ( i != 0 ) {
935  arg2 = AT.pWorkSpace[args];
936  AT.pWorkSpace[args] = arg1;
937  AT.pWorkSpace[args+1] = arg2;
938  }
939  m = mh = AT.WorkPointer;
940  mm = arg1+ARGHEAD; i = *mm;
941  NCOPY(m,mm,i);
942  AT.WorkPointer = m;
943  istart = 1;
944  argsdone++;
945  goto oneterm;
946  }
947  }
948 /*
949  #] There is a single term argument :
950  #[ Expand $ and expr :
951 
952  We have: 1: regular multiterm arguments
953  2: dollars
954  3: expressions.
955  The sum of them is numargs. Their addresses are in args. The problem is
956  that expansion will lead to allocations that we have to return and all
957  these allocations are different in nature.
958 */
959  action = 1;
960  abuf = (ARGBUFFER *)Malloc1(numargs*sizeof(ARGBUFFER),"argbuffer");
961  for ( i = 0; i < numargs; i++ ) {
962  arg1 = AT.pWorkSpace[args+i];
963  if ( *arg1 > 0 ) {
964  m = (WORD *)Malloc1(*arg1*sizeof(WORD),"argbuffer type 0");
965  abuf[i].buffer = m;
966  abuf[i].type = 0;
967  mm = arg1+ARGHEAD;
968  j = *arg1-ARGHEAD;
969  abuf[i].size = j;
970  if ( j ) argsdone++;
971  NCOPY(m,mm,j);
972  *m = 0;
973  }
974  else if ( *arg1 == -DOLLAREXPRESSION ) {
975  d = DolToTerms(BHEAD arg1[1]);
976  abuf[i].buffer = d->where;
977  abuf[i].type = 1;
978  abuf[i].dollar = d;
979  m = abuf[i].buffer;
980  if ( *m ) argsdone++;
981  while ( *m ) m+= *m;
982  abuf[i].size = m-abuf[i].buffer;
983  }
984  else if ( *arg1 == -EXPRESSION ) {
985  abuf[i].buffer = CreateExpression(BHEAD arg1[1]);
986  abuf[i].type = 2;
987  m = abuf[i].buffer;
988  if ( *m ) argsdone++;
989  while ( *m ) m+= *m;
990  abuf[i].size = m-abuf[i].buffer;
991  }
992  else {
993  MLOCK(ErrorMessageLock);
994  MesPrint("What argument is this?");
995  MUNLOCK(ErrorMessageLock);
996  goto CalledFrom;
997  }
998  }
999  for ( i = 0; i < numargs; i++ ) {
1000  arg1 = abuf[i].buffer;
1001  if ( *arg1 == 0 ) {}
1002  else if ( arg1[*arg1] == 0 ) {
1003 /*
1004  After expansion there is an argument with a single term
1005 */
1006  ab = abuf[i]; abuf[i] = abuf[0]; abuf[0] = ab;
1007  mh = abuf[0].buffer;
1008  for ( j = 1; j < numargs; j++ ) {
1009  m = abuf[j].buffer;
1010  while ( *m ) {
1011  GCDterms(BHEAD mh,m,mh); m += *m;
1012  argsdone++;
1013  if ( mh[0] == 4 && mh[1] == 1 && mh[2] == 1 && mh[3] == 3 ) {
1014  gcdisone = 1; sign = 1; break;
1015  }
1016  }
1017  if ( *m ) break;
1018  }
1019  mm = mh + *mh; if ( mm[-1] < 0 ) { sign = -1; mm[-1] = -mm[-1]; }
1020  mstop = mm - mm[-1]; m = mh+1; mlength = mm[-1];
1021  while ( tin < t ) *tout++ = *tin++;
1022  while ( m < mstop ) *tout++ = *m++;
1023  tin += tin[1];
1024  while ( tin < tstop ) *tout++ = *tin++;
1025  tlength = REDLENG(tlength);
1026  mlength = REDLENG(mlength);
1027  if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)mstop,mlength,
1028  (UWORD *)tout,&newlength) < 0 ) goto CalledFrom;
1029  mlength = INCLENG(newlength);
1030  tout += ABS(mlength);
1031  tout[-1] = mlength*sign;
1032  *termout = tout - termout;
1033  AT.WorkPointer = tout;
1034  if ( argsdone && Generator(BHEAD termout,level) < 0 ) goto CalledFrom;
1035  goto cleanup;
1036  }
1037  }
1038 /*
1039  There are only arguments with more than one term.
1040  We order them by size to make the computations as easy as possible.
1041 */
1042  for ( i = 1; i < numargs; i++ ) {
1043  for ( ii = i; ii > 0; ii-- ) {
1044  if ( abuf[ii-1].size <= abuf[ii].size ) break;
1045  ab = abuf[ii-1]; abuf[ii-1] = abuf[ii]; abuf[ii] = ab;
1046  }
1047  }
1048 /*
1049  #] Expand $ and expr :
1050  #[ Multiterm subexpressions :
1051 */
1052  ii = 0;
1053  gcdout = abuf[ii].buffer;
1054  for ( i = 0; i < numargs; i++ ) {
1055  if ( abuf[i].buffer[0] ) { gcdout = abuf[i].buffer; ii = i; i++; argsdone++; break; }
1056  }
1057  for ( ; i < numargs; i++ ) {
1058  if ( abuf[i].buffer[0] ) {
1059  g = GCDfunction3(BHEAD gcdout,abuf[i].buffer);
1060  argsdone++;
1061  if ( gcdout != abuf[ii].buffer ) M_free(gcdout,"gcdout");
1062  gcdout = g;
1063  if ( gcdout[*gcdout] == 0 && gcdout[0] == 4 && gcdout[1] == 1
1064  && gcdout[2] == 1 && gcdout[3] == 3 ) break;
1065  }
1066  }
1067  mm = gcdout;
1068 multiterms:;
1069  tlength = REDLENG(tlength);
1070  while ( *mm ) {
1071  tin = term; tout = termout; while ( tin < t ) *tout++ = *tin++;
1072  tin += t[1];
1073  mnext = mm + *mm; mlength = mnext[-1]; mstop = mnext - ABS(mlength);
1074  mm++;
1075  while ( mm < mstop ) *tout++ = *mm++;
1076  while ( tin < tstop ) *tout++ = *tin++;
1077  mlength = REDLENG(mlength);
1078  if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)mm,mlength,
1079  (UWORD *)tout,&newlength) < 0 ) goto CalledFrom;
1080  mlength = INCLENG(newlength);
1081  tout += ABS(mlength);
1082  tout[-1] = mlength;
1083  *termout = tout - termout;
1084  AT.WorkPointer = tout;
1085  if ( argsdone && Generator(BHEAD termout,level) < 0 ) goto CalledFrom;
1086  mm = mnext; /* next term */
1087  }
1088  if ( action && ( gcdout != abuf[ii].buffer ) ) M_free(gcdout,"gcdout");
1089 /*
1090  #] Multiterm subexpressions :
1091  #[ Cleanup :
1092 */
1093 cleanup:;
1094  if ( action ) {
1095  for ( i = 0; i < numargs; i++ ) {
1096  if ( abuf[i].type == 0 ) { M_free(abuf[i].buffer,"argbuffer type 0"); }
1097  else if ( abuf[i].type == 1 ) {
1098  d = abuf[i].dollar;
1099  if ( d->factors ) M_free(d->factors,"Dollar factors");
1100  M_free(d,"Copy of dollar variable");
1101  }
1102  else if ( abuf[i].type == 2 ) { M_free(abuf[i].buffer,"CreateExpression"); }
1103  }
1104  M_free(abuf,"argbuffer");
1105  }
1106 /*
1107  #] Cleanup :
1108 */
1109  AT.pWorkPointer = args;
1110  AT.WorkPointer = termout;
1111  return(0);
1112 
1113 CalledFrom:
1114  MLOCK(ErrorMessageLock);
1115  MesCall("GCDfunction");
1116  MUNLOCK(ErrorMessageLock);
1117  SETERROR(-1)
1118  return(-1);
1119 }
1120 
1121 /*
1122  #] GCDfunction :
1123  #[ GCDfunction3 :
1124 
1125  Finds the GCD of the two arguments which are buffers with terms.
1126  In principle the first buffer can have only one term.
1127 
1128  If both buffers have more than one term, we need to replace all
1129  non-symbolic objects by generated symbols and substitute that back
1130  afterwards. The rest we leave to the powerful routines.
1131  Philosophical problem: What do we do with GCD_(x/z+y,x+y*z) ?
1132 
1133  Method:
1134  If we have only negative powers of z and no positive powers we let
1135  the EXTRASYMBOLS do their job. When mixed, multiply the arguments with
1136  the negative powers with enough powers of z to eliminate the negative powers.
1137  The DENOMINATOR function is always eliminated with the mechanism as we
1138  cannot tell whether there are positive powers of its contents.
1139 */
1140 
1141 WORD *GCDfunction3(PHEAD WORD *in1, WORD *in2)
1142 {
1143  GETBIDENTITY
1144  WORD oldsorttype = AR.SortType, *ow = AT.WorkPointer;;
1145  WORD *t, *tt, *gcdout, *term1, *term2, *confree1, *confree2, *gcdout1, *proper1, *proper2;
1146  int i, actionflag1, actionflag2;
1147  WORD startebuf = cbuf[AT.ebufnum].numrhs;
1148  WORD tryterm1, tryterm2;
1149  if ( in2[*in2] == 0 ) { t = in1; in1 = in2; in2 = t; }
1150  if ( in1[*in1] == 0 ) { /* First input with only one term */
1151  gcdout = (WORD *)Malloc1((*in1+1)*sizeof(WORD),"gcdout");
1152  i = *in1; t = gcdout; tt = in1; NCOPY(t,tt,i); *t = 0;
1153  t = in2;
1154  while ( *t ) {
1155  GCDterms(BHEAD gcdout,t,gcdout);
1156  if ( gcdout[0] == 4 && gcdout[1] == 1
1157  && gcdout[2] == 1 && gcdout[3] == 3 ) break;
1158  t += *t;
1159  }
1160  AT.WorkPointer = ow;
1161  return(gcdout);
1162  }
1163 /*
1164  We need to take out the content from the two expressions
1165  and determine their GCD. This plays with the negative powers!
1166 */
1167  AR.SortType = SORTHIGHFIRST;
1168  term1 = TermMalloc("GCDfunction3-a");
1169  term2 = TermMalloc("GCDfunction3-b");
1170 
1171  confree1 = TakeContent(BHEAD in1,term1);
1172  tryterm1 = AN.tryterm; AN.tryterm = 0;
1173  confree2 = TakeContent(BHEAD in2,term2);
1174  tryterm2 = AN.tryterm; AN.tryterm = 0;
1175 /*
1176  confree1 = TakeSymbolContent(BHEAD in1,term1);
1177  confree2 = TakeSymbolContent(BHEAD in2,term2);
1178 */
1179  GCDterms(BHEAD term1,term2,term1);
1180  TermFree(term2,"GCDfunction3-b");
1181 /*
1182  Now we have to replace all non-symbols and symbols to a negative power
1183  by extra symbols.
1184 */
1185  if ( ( proper1 = PutExtraSymbols(BHEAD confree1,startebuf,&actionflag1) ) == 0 ) goto CalledFrom;
1186  if ( confree1 != in1 ) {
1187  if ( tryterm1 ) { TermFree(confree1,"TakeContent"); }
1188  else { M_free(confree1,"TakeContent"); }
1189  }
1190 /*
1191  TermFree(confree1,"TakeSymbolContent");
1192 */
1193  if ( ( proper2 = PutExtraSymbols(BHEAD confree2,startebuf,&actionflag2) ) == 0 ) goto CalledFrom;
1194  if ( confree2 != in2 ) {
1195  if ( tryterm2 ) { TermFree(confree2,"TakeContent"); }
1196  else { M_free(confree2,"TakeContent"); }
1197  }
1198 /*
1199  TermFree(confree2,"TakeSymbolContent");
1200 */
1201 /*
1202  And now the real work:
1203 */
1204  gcdout1 = poly_gcd(BHEAD proper1,proper2,0);
1205  M_free(proper1,"PutExtraSymbols");
1206  M_free(proper2,"PutExtraSymbols");
1207 
1208  AR.SortType = oldsorttype;
1209  if ( actionflag1 || actionflag2 ) {
1210  if ( ( gcdout = TakeExtraSymbols(BHEAD gcdout1,startebuf) ) == 0 ) goto CalledFrom;
1211  M_free(gcdout1,"gcdout");
1212  }
1213  else {
1214  gcdout = gcdout1;
1215  }
1216 
1217  cbuf[AT.ebufnum].numrhs = startebuf;
1218 /*
1219  Now multiply gcdout by term1
1220 */
1221  if ( term1[0] != 4 || term1[3] != 3 || term1[1] != 1 || term1[2] != 1 ) {
1222  AN.tryterm = -1;
1223  if ( ( gcdout1 = MultiplyWithTerm(BHEAD gcdout,term1,2) ) == 0 ) goto CalledFrom;
1224  AN.tryterm = 0;
1225  M_free(gcdout,"gcdout");
1226  gcdout = gcdout1;
1227  }
1228  TermFree(term1,"GCDfunction3-a");
1229  AT.WorkPointer = ow;
1230  return(gcdout);
1231 CalledFrom:
1232  AN.tryterm = 0;
1233  MLOCK(ErrorMessageLock);
1234  MesCall("GCDfunction3");
1235  MUNLOCK(ErrorMessageLock);
1236  return(0);
1237 }
1238 
1239 /*
1240  #] GCDfunction3 :
1241  #[ PutExtraSymbols :
1242 */
1243 
1244 WORD *PutExtraSymbols(PHEAD WORD *in,WORD startebuf,int *actionflag)
1245 {
1246  WORD *termout = AT.WorkPointer;
1247  int action;
1248  *actionflag = 0;
1249  NewSort(BHEAD0);
1250  while ( *in ) {
1251  if ( ( action = LocalConvertToPoly(BHEAD in,termout,startebuf,0) ) < 0 ) {
1252  LowerSortLevel();
1253  goto CalledFrom;
1254  }
1255  if ( action > 0 ) *actionflag = 1;
1256  StoreTerm(BHEAD termout);
1257  in += *in;
1258  }
1259  if ( EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 ) goto CalledFrom;
1260  return(termout);
1261 CalledFrom:
1262  MLOCK(ErrorMessageLock);
1263  MesCall("PutExtraSymbols");
1264  MUNLOCK(ErrorMessageLock);
1265  return(0);
1266 }
1267 
1268 /*
1269  #] PutExtraSymbols :
1270  #[ TakeExtraSymbols :
1271 */
1272 
1273 WORD *TakeExtraSymbols(PHEAD WORD *in,WORD startebuf)
1274 {
1275  CBUF *C = cbuf+AC.cbufnum;
1276  CBUF *CC = cbuf+AT.ebufnum;
1277  WORD *oldworkpointer = AT.WorkPointer, *termout;
1278 
1279  termout = AT.WorkPointer;
1280  NewSort(BHEAD0);
1281  while ( *in ) {
1282  if ( ConvertFromPoly(BHEAD in,termout,numxsymbol,CC->numrhs-startebuf+numxsymbol,startebuf-numxsymbol,1) <= 0 ) {
1283  LowerSortLevel();
1284  goto CalledFrom;
1285  }
1286  in += *in;
1287  AT.WorkPointer = termout + *termout;
1288 /*
1289  ConvertFromPoly leaves terms with subexpressions. Hence:
1290 */
1291  if ( Generator(BHEAD termout,C->numlhs) ) {
1292  LowerSortLevel();
1293  goto CalledFrom;
1294  }
1295  }
1296  AT.WorkPointer = oldworkpointer;
1297  if ( EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 ) goto CalledFrom;
1298  return(termout);
1299 
1300 CalledFrom:
1301  MLOCK(ErrorMessageLock);
1302  MesCall("TakeExtraSymbols");
1303  MUNLOCK(ErrorMessageLock);
1304  return(0);
1305 }
1306 
1307 /*
1308  #] TakeExtraSymbols :
1309  #[ MultiplyWithTerm :
1310 */
1311 
1312 WORD *MultiplyWithTerm(PHEAD WORD *in, WORD *term, WORD par)
1313 {
1314  WORD *termout, *t, *tt, *tstop, *ttstop;
1315  WORD length, length1, length2;
1316  WORD oldsorttype = AR.SortType;
1317  COMPARE oldcompareroutine = AR.CompareRoutine;
1318  AR.CompareRoutine = &CompareSymbols;
1319 
1320  if ( par == 0 || par == 2 ) AR.SortType = SORTHIGHFIRST;
1321  else AR.SortType = SORTLOWFIRST;
1322  termout = AT.WorkPointer;
1323  NewSort(BHEAD0);
1324  while ( *in ) {
1325  tt = termout + 1;
1326  tstop = in + *in; tstop -= ABS(tstop[-1]); t = in + 1;
1327  while ( t < tstop ) *tt++ = *t++;
1328  ttstop = term + *term; ttstop -= ABS(ttstop[-1]); t = term + 1;
1329  while ( t < ttstop ) *tt++ = *t++;
1330  length1 = REDLENG(in[*in-1]); length2 = REDLENG(term[*term-1]);
1331  if ( MulRat(BHEAD (UWORD *)tstop,length1,
1332  (UWORD *)ttstop,length2,(UWORD *)tt,&length) ) goto CalledFrom;
1333  length = INCLENG(length);
1334  tt += ABS(length); tt[-1] = length;
1335  *termout = tt - termout;
1336  SymbolNormalize(termout);
1337  StoreTerm(BHEAD termout);
1338  in += *in;
1339  }
1340  if ( par == 2 ) {
1341 /* if ( AN.tryterm == 0 ) AN.tryterm = 1; */
1342  AN.tryterm = 0; /* For now */
1343  if ( EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 ) goto CalledFrom;
1344  }
1345  else {
1346  if ( EndSort(BHEAD termout,1) < 0 ) goto CalledFrom;
1347  }
1348 
1349  AR.CompareRoutine = oldcompareroutine;
1350 
1351  AR.SortType = oldsorttype;
1352  return(termout);
1353 
1354 CalledFrom:
1355  MLOCK(ErrorMessageLock);
1356  MesCall("MultiplyWithTerm");
1357  MUNLOCK(ErrorMessageLock);
1358  return(0);
1359 }
1360 
1361 /*
1362  #] MultiplyWithTerm :
1363  #[ TakeContent :
1364 */
1376 WORD *TakeContent(PHEAD WORD *in, WORD *term)
1377 {
1378  GETBIDENTITY
1379  WORD *t, *tstop, *tcom, *tout, *tstore, *r, *rstop, *m, *mm, *w, *ww, *wterm;
1380  WORD *tnext, *tt, *tterm, code[2];
1381  WORD *inp, a, *den;
1382  int i, j, k, action = 0, sign;
1383  UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMbuffer2, *ap;
1384  WORD GCDlen, GCDlen2, LCMlen, LCMlen2, length, redlength, len1, len2;
1385  tout = tstore = term+1;
1386 /*
1387  #[ INDEX :
1388 */
1389  t = in;
1390  tnext = t + *t;
1391  tstop = tnext-ABS(tnext[-1]);
1392  t++;
1393  while ( t < tstop ) {
1394  if ( *t == INDEX ) {
1395  i = t[1]; NCOPY(tout,t,i); break;
1396  }
1397  else t += t[1];
1398  }
1399  if ( tout > tstore ) { /* There are indices in the first term */
1400  t = tnext;
1401  while ( *t ) {
1402  tnext = t + *t;
1403  rstop = tnext - ABS(tnext[-1]);
1404  r = t+1;
1405  if ( r == rstop ) goto noindices;
1406  while ( r < rstop ) {
1407  if ( *r != INDEX ) { r += r[1]; continue; }
1408  m = tstore+2;
1409  while ( m < tout ) {
1410  for ( i = 2; i < r[1]; i++ ) {
1411  if ( *m == r[i] ) break;
1412  if ( *m > r[i] ) continue;
1413  mm = m+1;
1414  while ( mm < tout ) { mm[-1] = mm[0]; mm++; }
1415  tout--; tstore[1]--; m--;
1416  break;
1417  }
1418  m++;
1419  }
1420  }
1421  if ( r >= rstop || tout <= tstore+2 ) {
1422  tout = tstore; break;
1423  }
1424  }
1425  if ( tout > tstore+2 ) { /* Now we have to take out what is in tstore */
1426  t = in; w = in;
1427  while ( *t ) {
1428  wterm = w;
1429  tnext = t + *t; t++; w++;
1430  while ( *t != INDEX ) { i = t[1]; NCOPY(w,t,i); }
1431  tt = t + t[1]; t += 2; r = tstore+2; ww = w; *w++ = INDEX; w++;
1432  while ( r < tout && t < tt ) {
1433  if ( *r > *t ) { *w++ = *t++; }
1434  else if ( *r == *t ) { r++; t++; }
1435  else goto CalledFrom;
1436  }
1437  if ( r < tout ) goto CalledFrom;
1438  while ( t < tt ) *w++ = *t++;
1439  ww[1] = w - ww;
1440  if ( ww[1] == 2 ) w = ww;
1441  while ( t < tnext ) *w++ = *t++;
1442  *wterm = w - wterm;
1443  }
1444  *w = 0;
1445  }
1446 noindices:
1447  tstore = tout;
1448  }
1449 /*
1450  #] INDEX :
1451  #[ VECTOR/DELTA :
1452 */
1453  code[0] = VECTOR; code[1] = DELTA;
1454  for ( k = 0; k < 2; k++ ) {
1455  t = in;
1456  tnext = t + *t;
1457  tstop = tnext-ABS(tnext[-1]);
1458  t++;
1459  while ( t < tstop ) {
1460  if ( *t == code[k] ) {
1461  i = t[1]; NCOPY(tout,t,i); break;
1462  }
1463  else t += t[1];
1464  }
1465  if ( tout > tstore ) { /* There are vectors in the first term */
1466  t = tnext;
1467  while ( *t ) {
1468  tnext = t + *t;
1469  rstop = tnext - ABS(tnext[-1]);
1470  r = t+1;
1471  if ( r == rstop ) { tstore = tout; goto novectors; }
1472  while ( r < rstop ) {
1473  if ( *r != code[k] ) { r += r[1]; continue; }
1474  m = tstore+2;
1475  while ( m < tout ) {
1476  for ( i = 2; i < r[1]; i += 2 ) {
1477  if ( *m == r[i] && m[1] == r[i+1] ) break;
1478  if ( *m > r[i] || ( *m == r[i] && m[1] > r[i+1] ) ) continue;
1479  mm = m+2;
1480  while ( mm < tout ) { mm[-2] = mm[0]; mm[-1] = mm[1]; mm += 2; }
1481  tout -= 2; tstore[1] -= 2; m -= 2;
1482  break;
1483  }
1484  m += 2;
1485  }
1486  }
1487  if ( r >= rstop || tout <= tstore+2 ) {
1488  tout = tstore; break;
1489  }
1490  }
1491  if ( tout > tstore+2 ) { /* Now we have to take out what is in tstore */
1492  t = in; w = in;
1493  while ( *t ) {
1494  wterm = w;
1495  tnext = t + *t; t++; w++;
1496  while ( *t != code[k] ) { i = t[1]; NCOPY(w,t,i); }
1497  tt = t + t[1]; t += 2; r = tstore+2; ww = w; *w++ = code[k]; w++;
1498  while ( r < tout && t < tt ) {
1499  if ( ( *r > *t ) || ( *r == *t && r[1] > t[1] ) )
1500  { *w++ = *t++; *w++ = *t++; }
1501  else if ( *r == *t && r[1] == t[1] ) { r += 2; t += 2; }
1502  else goto CalledFrom;
1503  }
1504  if ( r < tout ) goto CalledFrom;
1505  while ( t < tt ) *w++ = *t++;
1506  ww[1] = w - ww;
1507  if ( ww[1] == 2 ) w = ww;
1508  while ( t < tnext ) *w++ = *t++;
1509  *wterm = w - wterm;
1510  }
1511  *w = 0;
1512  }
1513  tstore = tout;
1514  }
1515  }
1516 novectors:;
1517 /*
1518  #] VECTOR/DELTA :
1519  #[ FUNCTIONS :
1520 */
1521  t = in;
1522  tnext = t + *t;
1523  tstop = tnext-ABS(tnext[-1]);
1524  t++;
1525  tcom = 0;
1526  while ( t < tstop ) {
1527  if ( *t >= FUNCTION ) {
1528  if ( functions[*t-FUNCTION].commute ) {
1529  if ( tcom == 0 ) { tcom = tstore; }
1530  else {
1531  for ( i = 0; i < t[1]; i++ ) {
1532  if ( t[i] != tcom[i] ) {
1533  MLOCK(ErrorMessageLock);
1534  MesPrint("GCD or factorization of more than one noncommuting object not allowed");
1535  MUNLOCK(ErrorMessageLock);
1536  goto CalledFrom;
1537  }
1538  }
1539  }
1540  }
1541  i = t[1]; NCOPY(tstore,t,i);
1542  }
1543  else t += t[1];
1544  }
1545  if ( tout > tstore ) {
1546  t = tnext;
1547  while ( *t ) {
1548  tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1549  if ( t == tstop ) goto nofunctions;
1550  r = tstore;
1551  while ( r < tout ) {
1552  tt = t;
1553  while ( tt < tstop ) {
1554  for ( i = 0; i < r[1]; i++ ) {
1555  if ( r[i] != tt[i] ) break;
1556  }
1557  if ( i == r[1] ) { r += r[1]; goto nextr1; }
1558  }
1559 /*
1560  Not encountered in this term. Scratch from list
1561 */
1562  m = r; mm = r + r[1];
1563  while ( mm < tout ) *m++ = *mm++;
1564  tout = m;
1565 nextr1:;
1566  }
1567  if ( tout <= tstore ) break;
1568  t += *t;
1569  }
1570  }
1571  if ( tout > tstore ) {
1572 /*
1573  Now we have one or more functions left that are common in all terms.
1574  Take them out. We do this one by one.
1575 */
1576  r = tstore;
1577  while ( r < tout ) {
1578  t = in; ww = in; w = ww+1;
1579  while ( *t ) {
1580  tnext = t + *t;
1581  t++;
1582  for(;;) {
1583  for ( i = 0; i < r[1]; i++ ) {
1584  if ( t[i] != r[i] ) {
1585  j = t[1]; NCOPY(w,t,j);
1586  break;
1587  }
1588  }
1589  if ( i == r[1] ) {
1590  t += t[1];
1591  while ( t < tnext ) *w++ = *t++;
1592  *ww = w - ww;
1593  break;
1594  }
1595  }
1596  }
1597  r += r[1];
1598  *w = 0;
1599  }
1600 nofunctions:
1601  tstore = tout;
1602  }
1603 /*
1604  #] FUNCTIONS :
1605  #[ SYMBOL :
1606 
1607  We make a list of symbols and their minimal powers.
1608  This includes negative powers. In the end we have to multiply by the
1609  inverse of this list. That takes out all negative powers and leaves
1610  things ready for further processing.
1611 */
1612  tterm = AT.WorkPointer; tt = tterm+1;
1613  tout[0] = SYMBOL; tout[1] = 2;
1614  t = in;
1615  tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1616  while ( t < tstop ) {
1617  if ( *t == SYMBOL ) {
1618  for ( i = 0; i < t[1]; i++ ) tout[i] = t[i];
1619  break;
1620  }
1621  t += t[1];
1622  }
1623  t = tnext;
1624  while ( *t ) {
1625  tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1626  if ( t == tstop ) {
1627  tout[1] = 2;
1628  break;
1629  }
1630  else {
1631  while ( t < tstop ) {
1632  if ( *t == SYMBOL ) {
1633  MergeSymbolLists(BHEAD tout,t,-1);
1634  break;
1635  }
1636  t += t[1];
1637  }
1638  t = tnext;
1639  }
1640  }
1641  if ( tout[1] > 2 ) {
1642  t = tout;
1643  tt[0] = t[0]; tt[1] = t[1];
1644  for ( i = 2; i < t[1]; i += 2 ) {
1645  tt[i] = t[i]; tt[i+1] = -t[i+1];
1646  }
1647  tt += tt[1];
1648  tout += tout[1];
1649  action++;
1650  }
1651 /*
1652  #] SYMBOL :
1653  #[ DOTPRODUCT :
1654 
1655  We make a list of dotproducts and their minimal powers.
1656  This includes negative powers. In the end we have to multiply by the
1657  inverse of this list. That takes out all negative powers and leaves
1658  things ready for further processing.
1659 */
1660  tout[0] = DOTPRODUCT; tout[1] = 2;
1661  t = in;
1662  while ( *t ) {
1663  tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
1664  if ( t == tstop ) {
1665  tout[1] = 2;
1666  break;
1667  }
1668  while ( t < tstop ) {
1669  if ( *t == DOTPRODUCT ) {
1670  MergeDotproductLists(BHEAD tout,t,-1);
1671  break;
1672  }
1673  t += t[1];
1674  }
1675  t = tnext;
1676  }
1677  if ( tout[1] > 2 ) {
1678  t = tout;
1679  tt[0] = t[0]; tt[1] = t[1];
1680  for ( i = 2; i < t[1]; i += 3 ) {
1681  tt[i] = t[i]; tt[i+1] = t[i+1]; tt[i+2] = -t[i+2];
1682  }
1683  tt += tt[1];
1684  tout += tout[1];
1685  action++;
1686  }
1687 /*
1688  #] DOTPRODUCT :
1689  #[ Coefficient :
1690 
1691  Now we have to collect the GCD of the numerators
1692  and the LCM of the denominators.
1693 */
1694  AT.WorkPointer = tt;
1695  if ( AN.cmod != 0 ) {
1696  WORD x, ix, ip;
1697  t = in; tnext = t + *t; tstop = tnext - ABS(tnext[-1]);
1698  x = tstop[0];
1699  if ( tnext[-1] < 0 ) x += AC.cmod[0];
1700  if ( GetModInverses(x,(WORD)(AN.cmod[0]),&ix,&ip) ) goto CalledFrom;
1701  *tout++ = x; *tout++ = 1; *tout++ = 3;
1702  *tt++ = ix; *tt++ = 1; *tt++ = 3;
1703  }
1704  else {
1705  GCDbuffer = NumberMalloc("MakeInteger");
1706  GCDbuffer2 = NumberMalloc("MakeInteger");
1707  LCMbuffer = NumberMalloc("MakeInteger");
1708  LCMbuffer2 = NumberMalloc("MakeInteger");
1709  t = in;
1710  tnext = t + *t; length = tnext[-1];
1711  if ( length < 0 ) { sign = -1; length = -length; }
1712  else { sign = 1; }
1713  tstop = tnext - length;
1714  redlength = (length-1)/2;
1715  for ( i = 0; i < redlength; i++ ) {
1716  GCDbuffer[i] = (UWORD)(tstop[i]);
1717  LCMbuffer[i] = (UWORD)(tstop[redlength+i]);
1718  }
1719  GCDlen = LCMlen = redlength;
1720  while ( GCDbuffer[GCDlen-1] == 0 ) GCDlen--;
1721  while ( LCMbuffer[LCMlen-1] == 0 ) LCMlen--;
1722  t = tnext;
1723  while ( *t ) {
1724  tnext = t + *t; length = ABS(tnext[-1]);
1725  tstop = tnext - length; redlength = (length-1)/2;
1726  len1 = len2 = redlength;
1727  den = tstop + redlength;
1728  while ( tstop[len1-1] == 0 ) len1--;
1729  while ( den[len2-1] == 0 ) len2--;
1730  if ( GCDlen == 1 && GCDbuffer[0] == 1 ) {}
1731  else {
1732  GcdLong(BHEAD (UWORD *)tstop,len1,GCDbuffer,GCDlen,GCDbuffer2,&GCDlen2);
1733  ap = GCDbuffer; GCDbuffer = GCDbuffer2; GCDbuffer2 = ap;
1734  a = GCDlen; GCDlen = GCDlen2; GCDlen2 = a;
1735  }
1736  if ( len2 == 1 && den[0] == 1 ) {}
1737  else {
1738  GcdLong(BHEAD LCMbuffer,LCMlen,(UWORD *)den,len2,LCMbuffer2,&LCMlen2);
1739  DivLong((UWORD *)den,len2,LCMbuffer2,LCMlen2,
1740  GCDbuffer2,&GCDlen2,(UWORD *)AT.WorkPointer,&a);
1741  MulLong(LCMbuffer,LCMlen,GCDbuffer2,GCDlen2,LCMbuffer2,&LCMlen2);
1742  ap = LCMbuffer; LCMbuffer = LCMbuffer2; LCMbuffer2 = ap;
1743  a = LCMlen; LCMlen = LCMlen2; LCMlen2 = a;
1744  }
1745  t = tnext;
1746  }
1747  if ( GCDlen != 1 || GCDbuffer[0] != 1 || LCMlen != 1 || LCMbuffer[0] != 1 ) {
1748  redlength = GCDlen; if ( LCMlen > GCDlen ) redlength = LCMlen;
1749  for ( i = 0; i < GCDlen; i++ ) *tout++ = (WORD)(GCDbuffer[i]);
1750  for ( ; i < redlength; i++ ) *tout++ = 0;
1751  for ( i = 0; i < LCMlen; i++ ) *tout++ = (WORD)(LCMbuffer[i]);
1752  for ( ; i < redlength; i++ ) *tout++ = 0;
1753  *tout++ = (2*redlength+1)*sign;
1754  for ( i = 0; i < LCMlen; i++ ) *tt++ = (WORD)(LCMbuffer[i]);
1755  for ( ; i < redlength; i++ ) *tt++ = 0;
1756  for ( i = 0; i < GCDlen; i++ ) *tt++ = (WORD)(GCDbuffer[i]);
1757  for ( ; i < redlength; i++ ) *tt++ = 0;
1758  *tt++ = (2*redlength+1)*sign;
1759  action++;
1760  }
1761  else {
1762  *tout++ = 1; *tout++ = 1; *tout++ = 3*sign;
1763  *tt++ = 1; *tt++ = 1; *tt++ = 3*sign;
1764  if ( sign != 1 ) action++;
1765  }
1766  *tout = 0;
1767  NumberFree(LCMbuffer2,"MakeInteger");
1768  NumberFree(LCMbuffer ,"MakeInteger");
1769  NumberFree(GCDbuffer2,"MakeInteger");
1770  NumberFree(GCDbuffer ,"MakeInteger");
1771  }
1772 /*
1773  #] Coefficient :
1774  #[ Multiply by the inverse content :
1775 */
1776  if ( action ) {
1777  *tterm = tt - tterm;
1778  AT.WorkPointer = tt;
1779  inp = MultiplyWithTerm(BHEAD in,tterm,2);
1780  AT.WorkPointer = tterm;
1781  in = inp;
1782  }
1783 /*
1784  #] Multiply by the inverse content :
1785 */
1786  *term = tout - term;
1787  AT.WorkPointer = tterm;
1788  return(in);
1789 CalledFrom:
1790  MLOCK(ErrorMessageLock);
1791  MesCall("TakeContent");
1792  MUNLOCK(ErrorMessageLock);
1793  return(0);
1794 }
1795 
1796 /*
1797  #] TakeContent :
1798  #[ MergeSymbolLists :
1799 
1800  Merges the extra list into the old.
1801  If par == -1 we take minimum powers
1802  If par == 1 we take maximum powers
1803  If par == 0 we take minimum of the absolute value of the powers
1804  if one is positive and the other negative we get zero.
1805  We assume that the symbols are in order in both lists
1806 */
1807 
1808 int MergeSymbolLists(PHEAD WORD *old, WORD *extra, int par)
1809 {
1810  GETBIDENTITY
1811  WORD *new = TermMalloc("MergeSymbolLists");
1812  WORD *t1, *t2, *fill;
1813  int i1,i2;
1814  fill = new + 2; *new = SYMBOL;
1815  i1 = old[1] - 2; i2 = extra[1] - 2;
1816  t1 = old + 2; t2 = extra + 2;
1817  switch ( par ) {
1818  case -1:
1819  while ( i1 > 0 && i2 > 0 ) {
1820  if ( *t1 > *t2 ) {
1821  if ( t2[1] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1822  else t2 += 2;
1823  i2 -= 2;
1824  }
1825  else if ( *t1 < *t2 ) {
1826  if ( t1[1] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1827  else t1 += 2;
1828  i1 -= 2;
1829  }
1830  else if ( t1[1] < t2[1] ) {
1831  *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1832  i1 -= 2; i2 -=2;
1833  }
1834  else {
1835  *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1836  i1 -= 2; i2 -=2;
1837  }
1838  }
1839  for ( ; i1 > 0; i1 -= 2 ) {
1840  if ( t1[1] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1841  else t1 += 2;
1842  }
1843  for ( ; i2 > 0; i2 -= 2 ) {
1844  if ( t2[1] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1845  else t2 += 2;
1846  }
1847  break;
1848  case 1:
1849  while ( i1 > 0 && i2 > 0 ) {
1850  if ( *t1 > *t2 ) {
1851  if ( t2[1] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1852  else t2 += 2;
1853  i2 -=2;
1854  }
1855  else if ( *t1 < *t2 ) {
1856  if ( t1[1] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1857  else t1 += 2;
1858  i1 -= 2;
1859  }
1860  else if ( t1[1] > t2[1] ) {
1861  *fill++ = *t1++; *fill++ = *t1++; t2 += 2;
1862  i1 -= 2; i2 -=2;
1863  }
1864  else {
1865  *fill++ = *t2++; *fill++ = *t2++; t1 += 2;
1866  i1 -= 2; i2 -=2;
1867  }
1868  }
1869  for ( ; i1 > 0; i1 -= 2 ) {
1870  if ( t1[1] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; }
1871  else t1 += 2;
1872  }
1873  for ( ; i2 > 0; i2 -= 2 ) {
1874  if ( t2[1] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; }
1875  else t2 += 2;
1876  }
1877  break;
1878  case 0:
1879  while ( i1 > 0 && i2 > 0 ) {
1880  if ( *t1 > *t2 ) {
1881  t2 += 2; i2 -= 2;
1882  }
1883  else if ( *t1 < *t2 ) {
1884  t1 += 2; i1 -= 2;
1885  }
1886  else if ( ( t1[1] > 0 ) && ( t2[1] < 0 ) ) { t1 += 2; t2 += 2; i1 -= 2; i2 -= 2; }
1887  else if ( ( t1[1] < 0 ) && ( t2[1] > 0 ) ) { t1 += 2; t2 += 2; i1 -= 2; i2 -= 2; }
1888  else if ( t1[1] > 0 ) {
1889  if ( t1[1] < t2[1] ) {
1890  *fill++ = *t1++; *fill++ = *t1++; t2 += 2; i2 -= 2;
1891  }
1892  else {
1893  *fill++ = *t2++; *fill++ = *t2++; t1 += 2; i1 -= 2;
1894  }
1895  }
1896  else {
1897  if ( t2[1] < t1[1] ) {
1898  *fill++ = *t2++; *fill++ = *t2++; t1 += 2; i1 -= 2; i2 -= 2;
1899  }
1900  else {
1901  *fill++ = *t1++; *fill++ = *t1++; t2 += 2; i1 -= 2; i2 -= 2;
1902  }
1903  }
1904  }
1905  for ( ; i1 > 0; i1-- ) *fill++ = *t1++;
1906  for ( ; i2 > 0; i2-- ) *fill++ = *t2++;
1907  break;
1908  }
1909  i1 = new[1] = fill - new;
1910  t2 = new; t1 = old; NCOPY(t1,t2,i1);
1911  TermFree(new,"MergeSymbolLists");
1912  return(0);
1913 }
1914 
1915 /*
1916  #] MergeSymbolLists :
1917  #[ MergeDotproductLists :
1918 
1919  Merges the extra list into the old.
1920  If par == -1 we take minimum powers
1921  If par == 1 we take maximum powers
1922  If par == 0 we take minimum of the absolute value of the powers
1923  if one is positive and the other negative we get zero.
1924  We assume that the dotproducts are in order in both lists
1925 */
1926 
1927 int MergeDotproductLists(PHEAD WORD *old, WORD *extra, int par)
1928 {
1929  GETBIDENTITY
1930  WORD *new = TermMalloc("MergeDotproductLists");
1931  WORD *t1, *t2, *fill;
1932  int i1,i2;
1933  fill = new + 2;
1934  i1 = old[1] - 2; i2 = extra[1] - 2;
1935  t1 = old + 2; t2 = extra + 2;
1936  switch ( par ) {
1937  case -1:
1938  while ( i1 > 0 && i2 > 0 ) {
1939  if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1940  if ( t2[2] < 0 ) { *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; }
1941  else t2 += 3;
1942  }
1943  else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1944  if ( t1[2] < 0 ) { *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; }
1945  else t1 += 3;
1946  }
1947  else if ( t1[2] < t2[2] ) {
1948  *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1949  }
1950  else {
1951  *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1952  }
1953  }
1954  break;
1955  case 1:
1956  while ( i1 > 0 && i2 > 0 ) {
1957  if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1958  if ( t2[2] > 0 ) { *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; }
1959  else t2 += 3;
1960  }
1961  else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1962  if ( t1[2] > 0 ) { *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; }
1963  else t1 += 3;
1964  }
1965  else if ( t1[2] > t2[2] ) {
1966  *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1967  }
1968  else {
1969  *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1970  }
1971  }
1972  break;
1973  case 0:
1974  while ( i1 > 0 && i2 > 0 ) {
1975  if ( ( *t1 > *t2 ) || ( *t1 == *t2 && t1[1] > t2[1] ) ) {
1976  t2 += 3;
1977  }
1978  else if ( ( *t1 < *t2 ) || ( *t1 == *t2 && t1[1] < t2[1] ) ) {
1979  t1 += 3;
1980  }
1981  else if ( ( t1[2] > 0 ) && ( t2[2] < 0 ) ) { t1 += 3; t2 += 3; }
1982  else if ( ( t1[2] < 0 ) && ( t2[2] > 0 ) ) { t1 += 3; t2 += 3; }
1983  else if ( t1[2] > 0 ) {
1984  if ( t1[2] < t2[2] ) {
1985  *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1986  }
1987  else {
1988  *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1989  }
1990  }
1991  else {
1992  if ( t2[2] < t1[2] ) {
1993  *fill++ = *t2++; *fill++ = *t2++; *fill++ = *t2++; t1 += 3;
1994  }
1995  else {
1996  *fill++ = *t1++; *fill++ = *t1++; *fill++ = *t1++; t2 += 3;
1997  }
1998  }
1999  }
2000  break;
2001  }
2002  i1 = new[1] = fill - new;
2003  t2 = new; t1 = old; NCOPY(t1,t2,i1);
2004  TermFree(new,"MergeDotproductLists");
2005  return(0);
2006 }
2007 
2008 /*
2009  #] MergeDotproductLists :
2010  #[ CreateExpression :
2011 
2012  Looks for the expression in the argument, reads it and puts it
2013  in a buffer. Returns the address of the buffer.
2014  We send the expression through the Generator system, because there
2015  may be unsubstituted (sub)expressions as in
2016  Local F = (a+b);
2017  Local G = gcd_(F,...);
2018 */
2019 
2020 WORD *CreateExpression(PHEAD WORD nexp)
2021 {
2022  GETBIDENTITY
2023  CBUF *C = cbuf+AC.cbufnum;
2024  POSITION startposition, oldposition;
2025  FILEHANDLE *fi;
2026  WORD *term, *oldipointer = AR.CompressPointer;
2027 ;
2028  switch ( Expressions[nexp].status ) {
2029  case HIDDENLEXPRESSION:
2030  case HIDDENGEXPRESSION:
2031  case DROPHLEXPRESSION:
2032  case DROPHGEXPRESSION:
2033  case UNHIDELEXPRESSION:
2034  case UNHIDEGEXPRESSION:
2035  AR.GetOneFile = 2; fi = AR.hidefile;
2036  break;
2037  default:
2038  AR.GetOneFile = 0; fi = AR.infile;
2039  break;
2040  }
2041  SeekScratch(fi,&oldposition);
2042  startposition = AS.OldOnFile[nexp];
2043  term = AT.WorkPointer;
2044  if ( GetOneTerm(BHEAD term,fi,&startposition,0) <= 0 ) goto CalledFrom;
2045  NewSort(BHEAD0);
2046  AR.CompressPointer = oldipointer;
2047  while ( GetOneTerm(BHEAD term,fi,&startposition,0) > 0 ) {
2048  AT.WorkPointer = term + *term;
2049  if ( Generator(BHEAD term,C->numlhs) ) {
2050  LowerSortLevel();
2051  goto CalledFrom;
2052  }
2053  AR.CompressPointer = oldipointer;
2054  }
2055  AT.WorkPointer = term;
2056  if ( EndSort(BHEAD (WORD *)((VOID *)(&term)),2) < 0 ) goto CalledFrom;
2057  SetScratch(fi,&oldposition);
2058  return(term);
2059 CalledFrom:
2060  MLOCK(ErrorMessageLock);
2061  MesCall("CreateExpression");
2062  MUNLOCK(ErrorMessageLock);
2063  Terminate(-1);
2064  return(0);
2065 }
2066 
2067 /*
2068  #] CreateExpression :
2069  #[ GCDterms : GCD of two terms
2070 
2071  Computes the GCD of two terms.
2072  Output in termout.
2073  termout may overlap with term1.
2074 */
2075 
2076 int GCDterms(PHEAD WORD *term1, WORD *term2, WORD *termout)
2077 {
2078  GETBIDENTITY
2079  WORD *t1, *t1stop, *t1next, *t2, *t2stop, *t2next, *tout, *tt1, *tt2;
2080  int count1, count2, i, ii, x1, sign;
2081  WORD length1, length2;
2082  t1 = term1 + *term1; t1stop = t1 - ABS(t1[-1]); t1 = term1+1;
2083  t2 = term2 + *term2; t2stop = t2 - ABS(t2[-1]); t2 = term2+1;
2084  tout = termout+1;
2085  while ( t1 < t1stop ) {
2086  t1next = t1 + t1[1];
2087  t2 = term2+1;
2088  if ( *t1 == SYMBOL ) {
2089  while ( t2 < t2stop && *t2 != SYMBOL ) t2 += t2[1];
2090  if ( t2 < t2stop && *t2 == SYMBOL ) {
2091  t2next = t2+t2[1];
2092  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2093  while ( tt1 < t1next && tt2 < t2next ) {
2094  if ( *tt1 < *tt2 ) tt1 += 2;
2095  else if ( *tt1 > *tt2 ) tt2 += 2;
2096  else if ( ( tt1[1] > 0 && tt2[1] < 0 ) ||
2097  ( tt2[1] > 0 && tt1[1] < 0 ) ) {
2098  tt1 += 2; tt2 += 2;
2099  }
2100  else {
2101  x1 = tt1[1];
2102  if ( tt1[1] < 0 ) { if ( tt2[1] > x1 ) x1 = tt2[1]; }
2103  else { if ( tt2[1] < x1 ) x1 = tt2[1]; }
2104  tout[count1+2] = *tt1;
2105  tout[count1+3] = x1;
2106  tt1 += 2; tt2 += 2;
2107  count1 += 2;
2108  }
2109  }
2110  if ( count1 > 0 ) {
2111  *tout = SYMBOL; tout[1] = count1+2; tout += tout[1];
2112  }
2113  }
2114  }
2115  else if ( *t1 == DOTPRODUCT ) {
2116  while ( t2 < t2stop && *t2 != DOTPRODUCT ) t2 += t2[1];
2117  if ( t2 < t2stop && *t2 == DOTPRODUCT ) {
2118  t2next = t2+t2[1];
2119  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2120  while ( tt1 < t1next && tt2 < t2next ) {
2121  if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 3;
2122  else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 3;
2123  else if ( ( tt1[2] > 0 && tt2[2] < 0 ) ||
2124  ( tt2[2] > 0 && tt1[2] < 0 ) ) {
2125  tt1 += 3; tt2 += 3;
2126  }
2127  else {
2128  x1 = tt1[2];
2129  if ( tt1[2] < 0 ) { if ( tt2[2] > x1 ) x1 = tt2[2]; }
2130  else { if ( tt2[2] < x1 ) x1 = tt2[2]; }
2131  tout[count1+2] = *tt1;
2132  tout[count1+3] = tt1[1];
2133  tout[count1+4] = x1;
2134  tt1 += 3; tt2 += 3;
2135  count1 += 3;
2136  }
2137  }
2138  if ( count1 > 0 ) {
2139  *tout = DOTPRODUCT; tout[1] = count1+2; tout += tout[1];
2140  }
2141  }
2142  }
2143  else if ( *t1 == VECTOR ) {
2144  while ( t2 < t2stop && *t2 != VECTOR ) t2 += t2[1];
2145  if ( t2 < t2stop && *t2 == VECTOR ) {
2146  t2next = t2+t2[1];
2147  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2148  while ( tt1 < t1next && tt2 < t2next ) {
2149  if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2150  else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2151  else {
2152  tout[count1+2] = *tt1;
2153  tout[count1+3] = tt1[1];
2154  tt1 += 2; tt2 += 2;
2155  count1 += 2;
2156  }
2157  }
2158  if ( count1 > 0 ) {
2159  *tout = VECTOR; tout[1] = count1+2; tout += tout[1];
2160  }
2161  }
2162  }
2163  else if ( *t1 == INDEX ) {
2164  while ( t2 < t2stop && *t2 != INDEX ) t2 += t2[1];
2165  if ( t2 < t2stop && *t2 == INDEX ) {
2166  t2next = t2+t2[1];
2167  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2168  while ( tt1 < t1next && tt2 < t2next ) {
2169  if ( *tt1 < *tt2 ) tt1 += 1;
2170  else if ( *tt1 > *tt2 ) tt2 += 1;
2171  else {
2172  tout[count1+2] = *tt1;
2173  tt1 += 1; tt2 += 1;
2174  count1 += 1;
2175  }
2176  }
2177  if ( count1 > 0 ) {
2178  *tout = INDEX; tout[1] = count1+2; tout += tout[1];
2179  }
2180  }
2181  }
2182  else if ( *t1 == DELTA ) {
2183  while ( t2 < t2stop && *t2 != DELTA ) t2 += t2[1];
2184  if ( t2 < t2stop && *t2 == DELTA ) {
2185  t2next = t2+t2[1];
2186  tt1 = t1+2; tt2 = t2+2; count1 = 0;
2187  while ( tt1 < t1next && tt2 < t2next ) {
2188  if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2189  else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2190  else {
2191  tout[count1+2] = *tt1;
2192  tout[count1+3] = tt1[1];
2193  tt1 += 2; tt2 += 2;
2194  count1 += 2;
2195  }
2196  }
2197  if ( count1 > 0 ) {
2198  *tout = DELTA; tout[1] = count1+2; tout += tout[1];
2199  }
2200  }
2201  }
2202  else if ( *t1 >= FUNCTION ) { /* noncommuting functions? Forbidden! */
2203 /*
2204  Count how many times this function occurs.
2205  Then count how many times it is in term2.
2206 */
2207  count1 = 1;
2208  while ( t1next < t1stop && *t1 == *t1next && t1[1] == t1next[1] ) {
2209  for ( i = 2; i < t1[1]; i++ ) {
2210  if ( t1[i] != t1next[i] ) break;
2211  }
2212  if ( i < t1[1] ) break;
2213  count1++;
2214  t1next += t1next[1];
2215  }
2216  count2 = 0;
2217  while ( t2 < t2stop ) {
2218  if ( *t2 == *t1 && t2[1] == t1[1] ) {
2219  for ( i = 2; i < t1[1]; i++ ) {
2220  if ( t2[i] != t1[i] ) break;
2221  }
2222  if ( i >= t1[1] ) count2++;
2223  }
2224  t2 += t2[1];
2225  }
2226  if ( count1 < count2 ) count2 = count1; /* number of common occurrences */
2227  if ( count2 > 0 ) {
2228  if ( tout == t1 ) {
2229  while ( count2 > 0 ) { tout += tout[1]; count2--; }
2230  }
2231  else {
2232  i = t1[1]*count2;
2233  NCOPY(tout,t1,i);
2234  }
2235  }
2236  }
2237  t1 = t1next;
2238  }
2239 /*
2240  Now the coefficients. They are in t1stop and t2stop. Should go to tout.
2241 */
2242  sign = 1;
2243  length1 = term1[*term1-1]; ii = i = ABS(length1); t1 = t1stop;
2244  if ( t1 != tout ) { NCOPY(tout,t1,i); tout -= ii; }
2245  length2 = term2[*term2-1];
2246  if ( length1 < 0 && length2 < 0 ) sign = -1;
2247  if ( AccumGCD(BHEAD (UWORD *)tout,&length1,(UWORD *)t2stop,length2) ) {
2248  MLOCK(ErrorMessageLock);
2249  MesCall("GCDterms");
2250  MUNLOCK(ErrorMessageLock);
2251  SETERROR(-1)
2252  }
2253  if ( sign < 0 && length1 > 0 ) length1 = -length1;
2254  tout += ABS(length1); tout[-1] = length1;
2255  *termout = tout - termout; *tout = 0;
2256  return(0);
2257 }
2258 
2259 /*
2260  #] GCDterms :
2261  #[ ReadPolyRatFun :
2262 */
2263 
2264 int ReadPolyRatFun(PHEAD WORD *term)
2265 {
2266  WORD *oldworkpointer = AT.WorkPointer;
2267  int flag, i;
2268  WORD *t, *fun, *nextt, *num, *den, *t1, *t2, size, numsize, densize;
2269  WORD *term1, *term2, *confree1, *confree2, *gcd, *num1, *den1, move, *newnum, *newden;
2270  WORD *tstop, *m1, *m2;
2271  WORD oldsorttype = AR.SortType;
2272  COMPARE oldcompareroutine = AR.CompareRoutine;
2273  AR.SortType = SORTHIGHFIRST;
2274  AR.CompareRoutine = &CompareSymbols;
2275 
2276  tstop = term + *term; tstop -= ABS(tstop[-1]);
2277  if ( term + *term == AT.WorkPointer ) flag = 1;
2278  else flag = 0;
2279  t = term+1;
2280  while ( t < tstop ) {
2281  if ( *t != AR.PolyFun ) { t += t[1]; continue; }
2282  if ( ( t[2] & MUSTCLEANPRF ) == 0 ) { t += t[1]; continue; }
2283  fun = t;
2284  nextt = t + t[1];
2285  if ( fun[1] > FUNHEAD && fun[FUNHEAD] == -SNUMBER && fun[FUNHEAD+1] == 0 )
2286  { *term = 0; break; }
2287  if ( FromPolyRatFun(BHEAD fun, &num, &den) > 0 ) { t = nextt; continue; }
2288  if ( *num == ARGHEAD ) { *term = 0; break; }
2289 /*
2290  Now we have num and den. Both are in general argument notation,
2291  but can also be used as expressions as in num+ARGHEAD, den+ARGHEAD.
2292  We need the gcd. For this we have to take out the contents
2293  because PreGCD does not like contents.
2294 */
2295  term1 = TermMalloc("ReadPolyRatFun");
2296  term2 = TermMalloc("ReadPolyRatFun");
2297  confree1 = TakeSymbolContent(BHEAD num+ARGHEAD,term1);
2298  confree2 = TakeSymbolContent(BHEAD den+ARGHEAD,term2);
2299  GCDclean(BHEAD term1,term2);
2300 /* gcd = PreGCD(BHEAD confree1,confree2,1); */
2301  gcd = poly_gcd(BHEAD confree1,confree2,1);
2302  newnum = PolyDiv(BHEAD confree1,gcd,"ReadPolyRatFun");
2303  newden = PolyDiv(BHEAD confree2,gcd,"ReadPolyRatFun");
2304  TermFree(confree2,"ReadPolyRatFun");
2305  TermFree(confree1,"ReadPolyRatFun");
2306  num1 = MULfunc(BHEAD term1,newnum);
2307  den1 = MULfunc(BHEAD term2,newden);
2308  TermFree(newnum,"ReadPolyRatFun");
2309  TermFree(newden,"ReadPolyRatFun");
2310 /* M_free(gcd,"poly_gcd"); */
2311  TermFree(gcd,"poly_gcd");
2312  TermFree(term1,"ReadPolyRatFun");
2313  TermFree(term2,"ReadPolyRatFun");
2314 /*
2315  Now we can put the function back together.
2316  Notice that we cannot use ToFast, because there is no reservation
2317  for the header of the argument. Fortunately there are only two
2318  types of fast arguments.
2319 */
2320  if ( num1[0] == 4 && num1[4] == 0 && num1[2] == 1 && num1[1] > 0 ) {
2321  numsize = 2; num1[0] = -SNUMBER;
2322  if ( num1[3] < 0 ) num1[1] = -num1[1];
2323  }
2324  else if ( num1[0] == 8 && num1[8] == 0 && num1[7] == 3 && num1[6] == 1
2325  && num1[5] == 1 && num1[1] == SYMBOL && num1[4] == 1 ) {
2326  numsize = 2; num1[0] = -SYMBOL; num1[1] = num1[3];
2327  }
2328  else { m1 = num1; while ( *m1 ) m1 += *m1; numsize = (m1-num1)+ARGHEAD; }
2329  if ( den1[0] == 4 && den1[4] == 0 && den1[2] == 1 && den1[1] > 0 ) {
2330  densize = 2; den1[0] = -SNUMBER;
2331  if ( den1[3] < 0 ) den1[1] = -den1[1];
2332  }
2333  else if ( den1[0] == 8 && den1[8] == 0 && den1[7] == 3 && den1[6] == 1
2334  && den1[5] == 1 && den1[1] == SYMBOL && den1[4] == 1 ) {
2335  densize = 2; den1[0] = -SYMBOL; den1[1] = den1[3];
2336  }
2337  else { m2 = den1; while ( *m2 ) m2 += *m2; densize = (m2-den1)+ARGHEAD; }
2338  size = FUNHEAD+numsize+densize;
2339 
2340  if ( size > fun[1] ) {
2341  move = size - fun[1];
2342  t1 = term+*term; t2 = t1+move;
2343  while ( t1 > nextt ) *--t2 = *--t1;
2344  tstop += move; nextt += move;
2345  *term += move;
2346  }
2347  else if ( size < fun[1] ) {
2348  move = fun[1]-size;
2349  t2 = fun+size; t1 = nextt;
2350  tstop -= move; nextt -= move;
2351  t = term+*term;
2352  while ( t1 < t ) *t2++ = *t1++;
2353  *term -= move;
2354  }
2355  else { /* no need to move anything */ }
2356  fun[1] = size; fun[2] = 0;
2357  t2 = fun+FUNHEAD; t1 = num1;
2358  if ( *num1 < 0 ) { *t2++ = num1[0]; *t2++ = num1[1]; }
2359  else { *t2++ = numsize; *t2++ = 0; FILLARG(t2);
2360  i = numsize-ARGHEAD; NCOPY(t2,t1,i) }
2361  t1 = den1;
2362  if ( *den1 < 0 ) { *t2++ = den1[0]; *t2++ = den1[1]; }
2363  else { *t2++ = densize; *t2++ = 0; FILLARG(t2);
2364  i = densize-ARGHEAD; NCOPY(t2,t1,i) }
2365 
2366  TermFree(num1,"MULfunc");
2367  TermFree(den1,"MULfunc");
2368  t = nextt;
2369  }
2370 
2371  if ( flag ) AT.WorkPointer = term +*term;
2372  else AT.WorkPointer = oldworkpointer;
2373  AR.CompareRoutine = oldcompareroutine;
2374  AR.SortType = oldsorttype;
2375  return(0);
2376 }
2377 
2378 /*
2379  #] ReadPolyRatFun :
2380  #[ FromPolyRatFun :
2381 */
2382 
2383 int FromPolyRatFun(PHEAD WORD *fun, WORD **numout, WORD **denout)
2384 {
2385  WORD *nextfun, *tt, *num, *den;
2386  int i;
2387  nextfun = fun + fun[1];
2388  fun += FUNHEAD;
2389  num = AT.WorkPointer;
2390  if ( *fun < 0 ) {
2391  if ( *fun != -SNUMBER && *fun != -SYMBOL ) goto Improper;
2392  ToGeneral(fun,num,0);
2393  tt = num + *num; *tt++ = 0;
2394  fun += 2;
2395  }
2396  else { i = *fun; tt = num; NCOPY(tt,fun,i); *tt++ = 0; }
2397  den = tt;
2398  if ( *fun < 0 ) {
2399  if ( *fun != -SNUMBER && *fun != -SYMBOL ) goto Improper;
2400  ToGeneral(fun,den,0);
2401  tt = den + *den; *tt++ = 0;
2402  fun += 2;
2403  }
2404  else { i = *fun; tt = den; NCOPY(tt,fun,i); *tt++ = 0; }
2405  *numout = num; *denout = den;
2406  if ( fun != nextfun ) { return(1); }
2407  AT.WorkPointer = tt;
2408  return(0);
2409 Improper:
2410  MLOCK(ErrorMessageLock);
2411  MesPrint("Improper use of PolyRatFun");
2412  MesCall("FromPolyRatFun");
2413  MUNLOCK(ErrorMessageLock);
2414  SETERROR(-1);
2415 }
2416 
2417 /*
2418  #] FromPolyRatFun :
2419  #[ TakeSymbolContent :
2420 */
2434 WORD *TakeSymbolContent(PHEAD WORD *in, WORD *term)
2435 {
2436  GETBIDENTITY
2437  WORD *t, *tstop, *tout, *tstore;
2438  WORD *tnext, *tt, *tterm;
2439  WORD *inp, a, *den, *oldworkpointer = AT.WorkPointer;
2440  int i, action = 0, sign, first;
2441  UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMbuffer2, *ap;
2442  WORD GCDlen, GCDlen2, LCMlen, LCMlen2, length, redlength, len1, len2;
2443  LONG j;
2444  tout = tstore = term+1;
2445 /*
2446  #[ SYMBOL :
2447 
2448  We make a list of symbols and their minimal powers.
2449  This includes negative powers. In the end we have to multiply by the
2450  inverse of this list. That takes out all negative powers and leaves
2451  things ready for further processing.
2452 */
2453  tterm = AT.WorkPointer; tt = tterm+1;
2454  tout[0] = SYMBOL; tout[1] = 2;
2455  t = in; first = 1;
2456  while ( *t ) {
2457  tnext = t + *t; tstop = tnext - ABS(tnext[-1]); t++;
2458  while ( t < tstop ) {
2459  if ( first ) {
2460  if ( *t == SYMBOL ) {
2461  for ( i = 0; i < t[1]; i++ ) tout[i] = t[i];
2462  goto didwork;
2463  }
2464  else {
2465  MLOCK(ErrorMessageLock);
2466  MesPrint ((char*)"ERROR: polynomials and polyratfuns must contain symbols only");
2467  MUNLOCK(ErrorMessageLock);
2468  Terminate(1);
2469  }
2470  }
2471  else if ( *t == SYMBOL ) {
2472  MergeSymbolLists(BHEAD tout,t,-1);
2473  goto didwork;
2474  }
2475  else {
2476  t += t[1];
2477  }
2478  }
2479 /*
2480  Here we come when there were no symbols. Only keep the negative ones.
2481 */
2482  if ( first == 0 ) {
2483  int j = 2;
2484  for ( i = 2; i < tout[1]; i += 2 ) {
2485  if ( tout[i+1] < 0 ) {
2486  if ( i == j ) { j += 2; }
2487  else { tout[j] = tout[i]; tout[j+1] = tout[i+1]; j += 2; }
2488  }
2489  }
2490  tout[1] = j;
2491  }
2492 didwork:;
2493  first = 0;
2494  t = tnext;
2495  }
2496  if ( tout[1] > 2 ) {
2497  t = tout;
2498  tt[0] = t[0]; tt[1] = t[1];
2499  for ( i = 2; i < t[1]; i += 2 ) {
2500  tt[i] = t[i]; tt[i+1] = -t[i+1];
2501  }
2502  tt += tt[1];
2503  tout += tout[1];
2504  action++;
2505  }
2506 /*
2507  #] SYMBOL :
2508  #[ Coefficient :
2509 
2510  Now we have to collect the GCD of the numerators
2511  and the LCM of the denominators.
2512 */
2513  AT.WorkPointer = tt;
2514  if ( AN.cmod != 0 ) {
2515  WORD x, ix, ip;
2516  t = in; tnext = t + *t; tstop = tnext - ABS(tnext[-1]);
2517  x = tstop[0];
2518  if ( tnext[-1] < 0 ) x += AC.cmod[0];
2519  if ( GetModInverses(x,(WORD)(AN.cmod[0]),&ix,&ip) ) goto CalledFrom;
2520  *tout++ = x; *tout++ = 1; *tout++ = 3;
2521  *tt++ = ix; *tt++ = 1; *tt++ = 3;
2522  }
2523  else {
2524  GCDbuffer = NumberMalloc("MakeInteger");
2525  GCDbuffer2 = NumberMalloc("MakeInteger");
2526  LCMbuffer = NumberMalloc("MakeInteger");
2527  LCMbuffer2 = NumberMalloc("MakeInteger");
2528  t = in;
2529  tnext = t + *t; length = tnext[-1];
2530  if ( length < 0 ) { sign = -1; length = -length; }
2531  else { sign = 1; }
2532  tstop = tnext - length;
2533  redlength = (length-1)/2;
2534  for ( i = 0; i < redlength; i++ ) {
2535  GCDbuffer[i] = (UWORD)(tstop[i]);
2536  LCMbuffer[i] = (UWORD)(tstop[redlength+i]);
2537  }
2538  GCDlen = LCMlen = redlength;
2539  while ( GCDbuffer[GCDlen-1] == 0 ) GCDlen--;
2540  while ( LCMbuffer[LCMlen-1] == 0 ) LCMlen--;
2541  t = tnext;
2542  while ( *t ) {
2543  tnext = t + *t; length = ABS(tnext[-1]);
2544  tstop = tnext - length; redlength = (length-1)/2;
2545  len1 = len2 = redlength;
2546  den = tstop + redlength;
2547  while ( tstop[len1-1] == 0 ) len1--;
2548  while ( den[len2-1] == 0 ) len2--;
2549  if ( GCDlen == 1 && GCDbuffer[0] == 1 ) {}
2550  else {
2551  GcdLong(BHEAD (UWORD *)tstop,len1,GCDbuffer,GCDlen,GCDbuffer2,&GCDlen2);
2552  ap = GCDbuffer; GCDbuffer = GCDbuffer2; GCDbuffer2 = ap;
2553  a = GCDlen; GCDlen = GCDlen2; GCDlen2 = a;
2554  }
2555  if ( len2 == 1 && den[0] == 1 ) {}
2556  else {
2557  GcdLong(BHEAD LCMbuffer,LCMlen,(UWORD *)den,len2,LCMbuffer2,&LCMlen2);
2558  DivLong((UWORD *)den,len2,LCMbuffer2,LCMlen2,
2559  GCDbuffer2,&GCDlen2,(UWORD *)AT.WorkPointer,&a);
2560  MulLong(LCMbuffer,LCMlen,GCDbuffer2,GCDlen2,LCMbuffer2,&LCMlen2);
2561  ap = LCMbuffer; LCMbuffer = LCMbuffer2; LCMbuffer2 = ap;
2562  a = LCMlen; LCMlen = LCMlen2; LCMlen2 = a;
2563  }
2564  t = tnext;
2565  }
2566  if ( GCDlen != 1 || GCDbuffer[0] != 1 || LCMlen != 1 || LCMbuffer[0] != 1 ) {
2567  redlength = GCDlen; if ( LCMlen > GCDlen ) redlength = LCMlen;
2568  for ( i = 0; i < GCDlen; i++ ) *tout++ = (WORD)(GCDbuffer[i]);
2569  for ( ; i < redlength; i++ ) *tout++ = 0;
2570  for ( i = 0; i < LCMlen; i++ ) *tout++ = (WORD)(LCMbuffer[i]);
2571  for ( ; i < redlength; i++ ) *tout++ = 0;
2572  *tout++ = (2*redlength+1)*sign;
2573  for ( i = 0; i < LCMlen; i++ ) *tt++ = (WORD)(LCMbuffer[i]);
2574  for ( ; i < redlength; i++ ) *tt++ = 0;
2575  for ( i = 0; i < GCDlen; i++ ) *tt++ = (WORD)(GCDbuffer[i]);
2576  for ( ; i < redlength; i++ ) *tt++ = 0;
2577  *tt++ = (2*redlength+1)*sign;
2578  action++;
2579  }
2580  else {
2581  *tout++ = 1; *tout++ = 1; *tout++ = 3*sign;
2582  *tt++ = 1; *tt++ = 1; *tt++ = 3*sign;
2583  if ( sign != 1 ) action++;
2584  }
2585  NumberFree(LCMbuffer2,"MakeInteger");
2586  NumberFree(LCMbuffer ,"MakeInteger");
2587  NumberFree(GCDbuffer2,"MakeInteger");
2588  NumberFree(GCDbuffer ,"MakeInteger");
2589  }
2590 /*
2591  #] Coefficient :
2592  #[ Multiply by the inverse content :
2593 */
2594  if ( action ) {
2595  *term = tout - term; *tout = 0;
2596  *tterm = tt - tterm; *tt = 0;
2597  AT.WorkPointer = tt;
2598  inp = MultiplyWithTerm(BHEAD in,tterm,2);
2599  AT.WorkPointer = tterm;
2600  t = inp; while ( *t ) t += *t;
2601  j = (t-inp); t = inp;
2602  if ( j*sizeof(WORD) > (size_t)(AM.MaxTer) ) goto OverWork;
2603  in = tout = TermMalloc("TakeSymbolContent");
2604  NCOPY(tout,t,j); *tout = 0;
2605  if ( AN.tryterm > 0 ) { TermFree(inp,"MultiplyWithTerm"); AN.tryterm = 0; }
2606  else { M_free(inp,"MultiplyWithTerm"); }
2607  }
2608  else {
2609  t = in; while ( *t ) t += *t;
2610  j = (t-in); t = in;
2611  if ( j*sizeof(WORD) > (size_t)(AM.MaxTer) ) goto OverWork;
2612  in = tout = TermMalloc("TakeSymbolContent");
2613  NCOPY(tout,t,j); *tout = 0;
2614  term[0] = 4; term[1] = 1; term[2] = 1; term[3] = 3; term[4] = 0;
2615  }
2616 /*
2617  #] Multiply by the inverse content :
2618  AT.WorkPointer = tterm + *tterm;
2619 */
2620  AT.WorkPointer = oldworkpointer;
2621 
2622  return(in);
2623 OverWork:
2624  MLOCK(ErrorMessageLock);
2625  MesPrint("Term too complex. Maybe increasing MaxTermSize can help");
2626  MUNLOCK(ErrorMessageLock);
2627 CalledFrom:
2628  MLOCK(ErrorMessageLock);
2629  MesCall("TakeSymbolContent");
2630  MUNLOCK(ErrorMessageLock);
2631  Terminate(-1);
2632  return(0);
2633 }
2634 
2635 /*
2636  #] TakeSymbolContent :
2637  #[ GCDclean :
2638 
2639  Takes a numerator and a denominator that each consist of a
2640  single term with only a coefficient and symbols and makes them
2641  into a proper fraction. Output overwrites input.
2642 */
2643 
2644 void GCDclean(PHEAD WORD *num, WORD *den)
2645 {
2646  WORD *out1 = TermMalloc("GCDclean");
2647  WORD *out2 = TermMalloc("GCDclean");
2648  WORD *t1, *t2, *r1, *r2, *t1stop, *t2stop, csize1, csize2, csize3, pow, sign;
2649  int i;
2650  t1stop = num+*num; sign = ( t1stop[-1] < 0 ) ? -1 : 1;
2651  csize1 = ABS(t1stop[-1]); t1stop -= csize1;
2652  t2stop = den+*den; if ( t2stop[-1] < 0 ) sign = -sign;
2653  csize2 = ABS(t2stop[-1]); t2stop -= csize2;
2654  t1 = num+1; t2 = den+1;
2655  r1 = out1+3; r2 = out2+3;
2656  if ( t1 == t1stop ) {
2657  if ( t2 < t2stop ) {
2658  for ( i = 2; i < t2[1]; i += 2 ) {
2659  if ( t2[i+1] < 0 ) { *r1++ = t2[i]; *r1++ = -t2[i+1]; }
2660  else { *r2++ = t2[i]; *r2++ = t2[i+1]; }
2661  }
2662  }
2663  }
2664  else if ( t2 == t2stop ) {
2665  for ( i = 2; i < t1[1]; i += 2 ) {
2666  if ( t1[i+1] < 0 ) { *r2++ = t1[i]; *r2++ = -t1[i+1]; }
2667  else { *r1++ = t1[i]; *r1++ = t1[i+1]; }
2668  }
2669  }
2670  else {
2671  t1 += 2; t2 += 2;
2672  while ( t1 < t1stop && t2 < t2stop ) {
2673  if ( *t1 < *t2 ) {
2674  if ( t1[1] > 0 ) { *r1++ = *t1; *r1++ = t1[1]; t1 += 2; }
2675  else if ( t1[1] < 0 ) { *r2++ = *t1; *r2++ = -t1[1]; t1 += 2; }
2676  }
2677  else if ( *t1 > *t2 ) {
2678  if ( t2[1] > 0 ) { *r2++ = *t2; *r2++ = t2[1]; t2 += 2; }
2679  else if ( t2[1] < 0 ) { *r1++ = *t2; *r1++ = -t2[1]; t2 += 2; }
2680  }
2681  else {
2682  pow = t1[1]-t2[1];
2683  if ( pow > 0 ) { *r1++ = *t1; *r1++ = pow; }
2684  else if ( pow < 0 ) { *r2++ = *t1; *r2++ = -pow; }
2685  t1 += 2; t2 += 2;
2686  }
2687  }
2688  while ( t1 < t1stop ) {
2689  if ( t1[1] < 0 ) { *r2++ = *t1; *r2++ = -t1[1]; }
2690  else { *r1++ = *t1; *r1++ = t1[1]; }
2691  t1 += 2;
2692  }
2693  while ( t2 < t2stop ) {
2694  if ( t2[1] < 0 ) { *r1++ = *t2; *r1++ = -t2[1]; }
2695  else { *r2++ = *t2; *r2++ = t2[1]; }
2696  t2 += 2;
2697  }
2698  }
2699  if ( r1 > out1+3 ) { out1[1] = SYMBOL; out1[2] = r1 - out1 - 1; }
2700  else r1 = out1+1;
2701  if ( r2 > out2+3 ) { out2[1] = SYMBOL; out2[2] = r2 - out2 - 1; }
2702  else r2 = out2+1;
2703 /*
2704  Now the coefficients.
2705 */
2706  csize1 = REDLENG(csize1);
2707  csize2 = REDLENG(csize2);
2708  if ( DivRat(BHEAD (UWORD *)t1stop,csize1,(UWORD *)t2stop,csize2,(UWORD *)r1,&csize3) ) {
2709  MLOCK(ErrorMessageLock);
2710  MesCall("GCDclean");
2711  MUNLOCK(ErrorMessageLock);
2712  Terminate(-1);
2713  }
2714  UnPack((UWORD *)r1,csize3,&csize2,&csize1);
2715  t2 = r1+ABS(csize3);
2716  for ( i = 0; i < csize2; i++ ) r2[i] = t2[i];
2717  r2 += csize2; *r2++ = 1;
2718  for ( i = 1; i < csize2; i++ ) *r2++ = 0;
2719  csize2 = INCLENG(csize2); *r2++ = csize2; *out2 = r2-out2;
2720  r1 += ABS(csize1); *r1++ = 1;
2721  for ( i = 1; i < ABS(csize1); i++ ) *r1++ = 0;
2722  csize1 = INCLENG(csize1); *r1++ = csize1; *out1 = r1-out1;
2723 
2724  t1 = num; t2 = out1; i = *out1; NCOPY(t1,t2,i); *t1 = 0;
2725  if ( sign < 0 ) t1[-1] = -t1[-1];
2726  t1 = den; t2 = out2; i = *out2; NCOPY(t1,t2,i); *t1 = 0;
2727 
2728  TermFree(out2,"GCDclean");
2729  TermFree(out1,"GCDclean");
2730 }
2731 
2732 /*
2733  #] GCDclean :
2734  #[ PolyDiv :
2735 
2736  Special stub function for polynomials that should fit inside a term.
2737  We make sure that the space is allocated by TermMalloc.
2738  This makes things much easier on the calling routines.
2739 */
2740 
2741 WORD *PolyDiv(PHEAD WORD *a,WORD *b,char *text)
2742 {
2743 /*
2744  Probably the following would work now
2745 */
2746  DUMMYUSE(text);
2747  return(poly_div(BHEAD a,b,1));
2748 /*
2749  WORD *quo, *qq;
2750  WORD *x, *xx;
2751  LONG i;
2752  quo = poly_div(BHEAD a,b,1);
2753  x = TermMalloc(text);
2754  qq = quo; while ( *qq ) qq += *qq;
2755  i = (qq-quo+1);
2756  if ( i*sizeof(WORD) > (size_t)(AM.MaxTer) ) {
2757  DUMMYUSE(text);
2758  MLOCK(ErrorMessageLock);
2759  MesPrint("PolyDiv: Term too complex. Maybe increasing MaxTermSize can help");
2760  MUNLOCK(ErrorMessageLock);
2761  Terminate(-1);
2762  }
2763  xx = x; qq = quo;
2764  NCOPY(xx,qq,i)
2765  TermFree(quo,"poly_div");
2766  return(x);
2767 */
2768 }
2769 
2770 /*
2771  #] PolyDiv :
2772  #] GCDfunction :
2773  #[ DIVfunction :
2774 
2775  Input: a div_ function that has two arguments inside a term.
2776  Action: Calculates [arg1/arg2] using polynomial techniques if needed.
2777  Output: The output result is combined with the remainder of the term
2778  and sent to Generator for further processing.
2779  Note that the output can be just a number or many terms.
2780  In case par == 0 the output is [arg1/arg2]
2781  In case par == 1 the output is [arg1%arg2]
2782  In case par == 2 the output is [inverse of arg1 modulus arg2]
2783  In case par == 3 the output is [arg1*arg2]
2784 */
2785 
2786 WORD divrem[4] = { DIVFUNCTION, REMFUNCTION, INVERSEFUNCTION, MULFUNCTION };
2787 
2788 int DIVfunction(PHEAD WORD *term,WORD level,int par)
2789 {
2790  GETBIDENTITY
2791  WORD *t, *tt, *r, *arg1 = 0, *arg2 = 0, *arg3 = 0, *termout;
2792  WORD *tstop, *tend, *r3, *rr, *rstop, tlength, rlength, newlength;
2793  WORD *proper1, *proper2, *proper3 = 0, numdol = -1;
2794  int numargs = 0, type1, type2, actionflag1, actionflag2;
2795  WORD startebuf = cbuf[AT.ebufnum].numrhs;
2796  int division = ( par <= 2 ); /* false for mul_ */
2797  if ( par < 0 || par > 3 ) {
2798  MLOCK(ErrorMessageLock);
2799  MesPrint("Internal error. Illegal parameter %d in DIVfunction.",par);
2800  MUNLOCK(ErrorMessageLock);
2801  Terminate(-1);
2802  }
2803 /*
2804  Find the function
2805 */
2806  tend = term + *term; tstop = tend - ABS(tend[-1]);
2807  t = term+1;
2808  while ( t < tstop ) {
2809  if ( *t != divrem[par] ) { t += t[1]; continue; }
2810  r = t + FUNHEAD;
2811  tt = t + t[1]; numargs = 0;
2812  while ( r < tt ) {
2813  if ( numargs == 0 ) { arg1 = r; }
2814  if ( numargs == 1 ) { arg2 = r; }
2815  if ( numargs == 2 && *r == -DOLLAREXPRESSION ) { numdol = r[1]; }
2816  numargs++;
2817  NEXTARG(r);
2818  }
2819  if ( numargs == 2 ) break;
2820  if ( division && numargs == 3 ) break;
2821  t = tt;
2822  }
2823  if ( t >= tstop ) {
2824  MLOCK(ErrorMessageLock);
2825  MesPrint("Internal error. Indicated div_ or rem_ function not encountered.");
2826  MUNLOCK(ErrorMessageLock);
2827  Terminate(-1);
2828  }
2829 /*
2830  We have two arguments in arg1 and arg2.
2831 */
2832  if ( division && *arg1 == -SNUMBER && arg1[1] == 0 ) {
2833  if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2834 zerozero:;
2835  MLOCK(ErrorMessageLock);
2836  MesPrint("0/0 in either div_ or rem_ function.");
2837  MUNLOCK(ErrorMessageLock);
2838  Terminate(-1);
2839  }
2840  if ( numdol >= 0 ) PutTermInDollar(0,numdol);
2841  return(0);
2842  }
2843  if ( division && *arg2 == -SNUMBER && arg2[1] == 0 ) {
2844 divzero:;
2845  MLOCK(ErrorMessageLock);
2846  MesPrint("Division by zero in either div_ or rem_ function.");
2847  MUNLOCK(ErrorMessageLock);
2848  Terminate(-1);
2849  }
2850  if ( !division ) {
2851  if ( (*arg1 == -SNUMBER && arg1[1] == 0) ||
2852  (*arg2 == -SNUMBER && arg2[1] == 0) ) {
2853  return(0);
2854  }
2855  }
2856  if ( ( arg1 = ConvertArgument(BHEAD arg1, &type1) ) == 0 ) goto CalledFrom;
2857  if ( ( arg2 = ConvertArgument(BHEAD arg2, &type2) ) == 0 ) goto CalledFrom;
2858  if ( division && *arg1 == 0 ) {
2859  if ( *arg2 == 0 ) {
2860  M_free(arg2,"DIVfunction");
2861  M_free(arg1,"DIVfunction");
2862  goto zerozero;
2863  }
2864  M_free(arg2,"DIVfunction");
2865  M_free(arg1,"DIVfunction");
2866  if ( numdol >= 0 ) PutTermInDollar(0,numdol);
2867  return(0);
2868  }
2869  if ( division && *arg2 == 0 ) {
2870  M_free(arg2,"DIVfunction");
2871  M_free(arg1,"DIVfunction");
2872  goto divzero;
2873  }
2874  if ( !division && (*arg1 == 0 || *arg2 == 0) ) {
2875  M_free(arg2,"DIVfunction");
2876  M_free(arg1,"DIVfunction");
2877  return(0);
2878  }
2879  if ( ( proper1 = PutExtraSymbols(BHEAD arg1,startebuf,&actionflag1) ) == 0 ) goto CalledFrom;
2880  if ( ( proper2 = PutExtraSymbols(BHEAD arg2,startebuf,&actionflag2) ) == 0 ) goto CalledFrom;
2881 /*
2882  if ( type2 == 0 ) M_free(arg2,"DIVfunction");
2883  else {
2884  DOLLARS d = ((DOLLARS)arg2)-1;
2885  if ( d->factors ) M_free(d->factors,"Dollar factors");
2886  M_free(d,"Copy of dollar variable");
2887  }
2888 */
2889  M_free(arg2,"DIVfunction");
2890 /*
2891  if ( type1 == 0 ) M_free(arg1,"DIVfunction");
2892  else {
2893  DOLLARS d = ((DOLLARS)arg1)-1;
2894  if ( d->factors ) M_free(d->factors,"Dollar factors");
2895  M_free(d,"Copy of dollar variable");
2896  }
2897 */
2898  M_free(arg1,"DIVfunction");
2899  if ( par == 0 ) proper3 = poly_div(BHEAD proper1, proper2,0);
2900  else if ( par == 1 ) proper3 = poly_rem(BHEAD proper1, proper2,0);
2901  else if ( par == 2 ) proper3 = poly_inverse(BHEAD proper1, proper2);
2902  else if ( par == 3 ) proper3 = poly_mul(BHEAD proper1, proper2);
2903  if ( proper3 == 0 ) goto CalledFrom;
2904  if ( actionflag1 || actionflag2 ) {
2905  if ( ( arg3 = TakeExtraSymbols(BHEAD proper3,startebuf) ) == 0 ) goto CalledFrom;
2906  M_free(proper3,"DIVfunction");
2907  }
2908  else {
2909  arg3 = proper3;
2910  }
2911  M_free(proper2,"DIVfunction");
2912  M_free(proper1,"DIVfunction");
2913  cbuf[AT.ebufnum].numrhs = startebuf;
2914  if ( *arg3 ) {
2915  termout = AT.WorkPointer;
2916  tlength = tend[-1];
2917  tlength = REDLENG(tlength);
2918  r3 = arg3;
2919  while ( *r3 ) {
2920  tt = term + 1; rr = termout + 1;
2921  while ( tt < t ) *rr++ = *tt++;
2922  r = r3 + 1;
2923  r3 = r3 + *r3;
2924  rstop = r3 - ABS(r3[-1]);
2925  while ( r < rstop ) *rr++ = *r++;
2926  tt += t[1];
2927  while ( tt < tstop ) *rr++ = *tt++;
2928  rlength = r3[-1];
2929  rlength = REDLENG(rlength);
2930  if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)rstop,rlength,
2931  (UWORD *)rr,&newlength) < 0 ) goto CalledFrom;
2932  rlength = INCLENG(newlength);
2933  rr += ABS(rlength);
2934  rr[-1] = rlength;
2935  *termout = rr - termout;
2936  AT.WorkPointer = rr;
2937  if ( Generator(BHEAD termout,level) ) goto CalledFrom;
2938  }
2939  AT.WorkPointer = termout;
2940  }
2941  M_free(arg3,"DIVfunction");
2942  return(0);
2943 CalledFrom:
2944  MLOCK(ErrorMessageLock);
2945  MesCall("DIVfunction");
2946  MUNLOCK(ErrorMessageLock);
2947  SETERROR(-1)
2948 }
2949 
2950 /*
2951  #] DIVfunction :
2952  #[ MULfunc :
2953 
2954  Multiplies two polynomials and puts the results in TermMalloc space.
2955 */
2956 
2957 WORD *MULfunc(PHEAD WORD *p1, WORD *p2)
2958 {
2959  WORD *prod,size1,size2,size3,*t,*tfill,*ps1,*ps2,sign1,sign2, error, *p3;
2960  UWORD *num1, *num2, *num3;
2961  int i;
2962  WORD oldsorttype = AR.SortType;
2963  COMPARE oldcompareroutine = AR.CompareRoutine;
2964  AR.SortType = SORTHIGHFIRST;
2965  AR.CompareRoutine = &CompareSymbols;
2966  num3 = NumberMalloc("MULfunc");
2967  prod = TermMalloc("MULfunc");
2968  NewSort(BHEAD0);
2969  while ( *p1 ) {
2970  ps1 = p1+*p1; num1 = (UWORD *)(ps1 - ABS(ps1[-1])); size1 = ps1[-1];
2971  if ( size1 < 0 ) { sign1 = -1; size1 = -size1; }
2972  else sign1 = 1;
2973  size1 = (size1-1)/2;
2974  p3 = p2;
2975  while ( *p3 ) {
2976  ps2 = p3+*p3; num2 = (UWORD *)(ps2 - ABS(ps2[-1])); size2 = ps2[-1];
2977  if ( size2 < 0 ) { sign2 = -1; size2 = -size2; }
2978  else sign2 = 1;
2979  size2 = (size2-1)/2;
2980  if ( MulLong(num1,size1,num2,size2,num3,&size3) ) {
2981  error = 1;
2982 CalledFrom:
2983  MLOCK(ErrorMessageLock);
2984  MesPrint(" Error %d",error);
2985  MesCall("MulFunc");
2986  MUNLOCK(ErrorMessageLock);
2987  Terminate(-1);
2988  }
2989  tfill = prod+1;
2990  t = p1+1; while ( t < (WORD *)num1 ) *tfill++ = *t++;
2991  t = p3+1; while ( t < (WORD *)num2 ) *tfill++ = *t++;
2992  t = (WORD *)num3;
2993  for ( i = 0; i < size3; i++ ) *tfill++ = *t++;
2994  *tfill++ = 1;
2995  for ( i = 1; i < size3; i++ ) *tfill++ = 0;
2996  *tfill++ = (2*size3+1)*sign1*sign2;
2997  prod[0] = tfill - prod;
2998  if ( SymbolNormalize(prod) ) { error = 2; goto CalledFrom; }
2999  if ( StoreTerm(BHEAD prod) ) { error = 3; goto CalledFrom; }
3000  p3 += *p3;
3001  }
3002  p1 += *p1;
3003  }
3004  NumberFree(num3,"MULfunc");
3005  EndSort(BHEAD prod,1);
3006  AR.CompareRoutine = oldcompareroutine;
3007  AR.SortType = oldsorttype;
3008  return(prod);
3009 }
3010 
3011 /*
3012  #] MULfunc :
3013  #[ ConvertArgument :
3014 
3015  Converts an argument to a general notation in allocated space.
3016 */
3017 
3018 WORD *ConvertArgument(PHEAD WORD *arg, int *type)
3019 {
3020  WORD *output, *t, *r;
3021  int i, size;
3022  if ( *arg > 0 ) {
3023  output = (WORD *)Malloc1((*arg)*sizeof(WORD),"ConvertArgument");
3024  i = *arg - ARGHEAD; t = arg + ARGHEAD; r = output;
3025  NCOPY(r,t,i);
3026  *r = 0; *type = 0;
3027  return(output);
3028  }
3029  if ( *arg == -EXPRESSION ) {
3030  *type = 0;
3031  return(CreateExpression(BHEAD arg[1]));
3032  }
3033  if ( *arg == -DOLLAREXPRESSION ) {
3034  DOLLARS d;
3035  *type = 1;
3036  d = DolToTerms(BHEAD arg[1]);
3037 /*
3038  The problem is that DolToTerms creates a copy of the dollar variable.
3039  If we just return d->where we create a memory leak. Hence we have to
3040  copy the contents of d->where to a new buffer
3041 */
3042  output = (WORD *)Malloc1((d->size+1)*sizeof(WORD),"Copy of dollar content");
3043  WCOPY(output,d->where,d->size+1);
3044  if ( d->factors ) { M_free(d->factors,"Dollar factors"); d->factors = 0; }
3045  M_free(d,"Copy of dollar variable");
3046  return(output);
3047  }
3048 #if ( FUNHEAD > 4 )
3049  size = FUNHEAD+5;
3050 #else
3051  size = 9;
3052 #endif
3053  output = (WORD *)Malloc1(size*sizeof(WORD),"ConvertArgument");
3054  switch(*arg) {
3055  case -SYMBOL:
3056  output[0] = 8; output[1] = SYMBOL; output[2] = 4; output[3] = arg[1];
3057  output[4] = 1; output[5] = 1; output[6] = 1; output[7] = 3;
3058  output[8] = 0;
3059  break;
3060  case -INDEX:
3061  case -VECTOR:
3062  case -MINVECTOR:
3063  output[0] = 7; output[1] = INDEX; output[2] = 3; output[3] = arg[1];
3064  output[4] = 1; output[5] = 1;
3065  if ( *arg == -MINVECTOR ) output[6] = -3;
3066  else output[6] = 3;
3067  output[7] = 0;
3068  break;
3069  case -SNUMBER:
3070  output[0] = 4;
3071  if ( arg[1] < 0 ) {
3072  output[1] = -arg[1]; output[2] = 1; output[3] = -3;
3073  }
3074  else {
3075  output[1] = arg[1]; output[2] = 1; output[3] = 3;
3076  }
3077  output[4] = 0;
3078  break;
3079  default:
3080  output[0] = FUNHEAD+4;
3081  output[1] = -*arg;
3082  output[2] = FUNHEAD;
3083  for ( i = 3; i <= FUNHEAD; i++ ) output[i] = 0;
3084  output[FUNHEAD+1] = 1;
3085  output[FUNHEAD+2] = 1;
3086  output[FUNHEAD+3] = 3;
3087  output[FUNHEAD+4] = 0;
3088  break;
3089  }
3090  *type = 0;
3091  return(output);
3092 }
3093 
3094 /*
3095  #] ConvertArgument :
3096  #[ ExpandRat :
3097 
3098  Expands the denominator of a PolyRatFun in the variable PolyFunVar.
3099  The output is a polyratfun with a single argument.
3100  In the case that there is a polyratfun with more than one argument
3101  or the dirtyflag is on, the argument(s) is/are normalized.
3102  The output overwrites the input.
3103 */
3104 
3105 char *TheErrorMessage[] = {
3106  "PolyRatFun not of a type that FORM will expand: incorrect variable inside."
3107  ,"Division by zero in PolyRatFun encountered in ExpandRat."
3108  ,"Irregular code in PolyRatFun encountered in ExpandRat."
3109  ,"Called from ExpandRat."
3110  ,"WorkSpace overflow. Change parameter WorkSpace in setup file?"
3111  };
3112 
3113 int ExpandRat(PHEAD WORD *fun)
3114 {
3115  WORD *r, *rr, *rrr, *tt, *tnext, *arg1, *arg2, *rmin = 0, *rmininv;
3116  WORD *rcoef, rsize, rcopy, *ow = AT.WorkPointer;
3117  WORD *numerator, *denominator, *rnext;
3118  WORD *thecopy, *rc, ncoef, newcoef, *m, *mm, nco, *outarg = 0;
3119  UWORD co[2], co1[2], co2[2];
3120  WORD OldPolyFunPow = AR.PolyFunPow;
3121  int i, j, minpow = 0, eppow, first, error = 0, ipoly;
3122  if ( fun[1] == FUNHEAD ) { return(0); }
3123  tnext = fun + fun[1];
3124  if ( fun[1] == fun[FUNHEAD]+FUNHEAD ) { /* Single argument */
3125  if ( fun[2] == 0 ) { goto done; }
3126 /*
3127  We have to normalize the argument. This could make it shorter.
3128 */
3129 NormArg:;
3130  if ( outarg == 0 ) outarg = TermMalloc("ExpandRat")+ARGHEAD;
3131  AT.TrimPower = 1;
3132  NewSort(BHEAD0);
3133  r = fun+FUNHEAD+ARGHEAD;
3134  if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3135  WORD minpow2 = MAXPOWER, *rrm;
3136  rrm = r;
3137  while ( rrm < tnext ) {
3138  if ( *rrm == 4 ) {
3139  if ( minpow2 > 0 ) minpow2 = 0;
3140  }
3141  else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3142  if ( minpow2 > 0 ) minpow2 = 0;
3143  }
3144  else {
3145  if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3146  if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3147  }
3148  else {
3149  MesPrint("Illegal term in expanded polyratfun.");
3150  goto onerror;
3151  }
3152  }
3153  rrm += *rrm;
3154  }
3155  AR.PolyFunPow += minpow2;
3156  }
3157  while ( r < tnext ) {
3158  rr = r + *r;
3159  i = *r; rrr = outarg; NCOPY(rrr,r,i);
3160  Normalize(BHEAD outarg);
3161  if ( *outarg > 0 ) StoreTerm(BHEAD outarg);
3162  }
3163  r = fun+FUNHEAD+ARGHEAD;
3164  EndSort(BHEAD r,1);
3165  AT.TrimPower = 0;
3166  if ( *r == 0 ) {
3167  fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0;
3168  fun[1] = FUNHEAD+2;
3169  }
3170  else {
3171  rr = fun+FUNHEAD;
3172  if ( ToFast(rr,rr) ) {
3173  NEXTARG(rr); fun[1] = rr - fun;
3174  }
3175  else {
3176  while ( *r ) r += *r;
3177  *rr = r-rr; rr[1] = CLEANFLAG;
3178  fun[1] = r - fun;
3179  }
3180  }
3181  fun[2] = CLEANFLAG;
3182  goto done;
3183  }
3184 /*
3185  First test whether we have only AR.PolyFunVar in the denominator
3186 */
3187  tt = fun + FUNHEAD;
3188  arg1 = arg2 = 0;
3189  if ( tt < tnext ) {
3190  arg1 = tt; NEXTARG(tt);
3191  if ( tt < tnext ) {
3192  arg2 = tt; NEXTARG(tt);
3193  if ( tt != tnext ) { arg1 = arg2 = 0; } /* more than two arguments */
3194  }
3195  }
3196  if ( arg2 == 0 ) {
3197  if ( *arg1 < 0 ) { fun[2] = CLEANFLAG; goto done; }
3198  if ( fun[2] == CLEANFLAG ) goto done;
3199  goto NormArg; /* Note: should not come here */
3200  }
3201 /*
3202  Produce the output argument in outarg
3203 */
3204  if ( outarg == 0 ) outarg = TermMalloc("ExpandRat")+ARGHEAD;
3205 
3206  if ( *arg2 <= 0 ) {
3207 /*
3208  These cases are trivial.
3209  We try as much as possible to write the output directly into the
3210  function. We just have to be extremely careful not to overwrite
3211  relevant information before we are finished with it.
3212 */
3213  if ( *arg2 == -SYMBOL && arg2[1] == AR.PolyFunVar ) {
3214  rr = r = fun+FUNHEAD+ARGHEAD;
3215  if ( *arg1 < 0 ) {
3216  if ( *arg1 == -SYMBOL ) {
3217  if ( arg1[1] == AR.PolyFunVar ) {
3218  *r++ = 4; *r++ = 1; *r++ = 1; *r++ = 3; *r = 0;
3219  }
3220  else {
3221  *r++ = 10; *r++ = SYMBOL; *r++ = 6;
3222  *r++ = arg1[1]; *r++ = 1;
3223  *r++ = AR.PolyFunVar; *r++ = -1;
3224  *r++ = 1; *r++ = 1; *r++ = 3; *r = 0;
3225  Normalize(BHEAD rr);
3226  }
3227  }
3228  else if ( *arg1 == -SNUMBER ) {
3229  nco = arg1[1];
3230  if ( nco == 0 ) { *r++ = 0; }
3231  else {
3232  *r++ = 8; *r++ = SYMBOL; *r++ = 4;
3233  *r++ = AR.PolyFunVar; *r++ = -1;
3234  *r++ = ABS(nco); *r++ = 1;
3235  if ( nco < 0 ) *r++ = -3;
3236  else *r++ = 3;
3237  *r = 0;
3238  }
3239  }
3240  else { error = 2; goto onerror; } /* should not happen! */
3241  }
3242  else { /* Multi-term numerator. */
3243  m = arg1+ARGHEAD;
3244  NewSort(BHEAD0); /* Technically maybe not needed */
3245  if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3246  WORD minpow2 = MAXPOWER, *rrm;
3247  rrm = m;
3248  while ( rrm < arg2 ) {
3249  if ( *rrm == 4 ) {
3250  if ( minpow2 > 0 ) minpow2 = 0;
3251  }
3252  else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3253  if ( minpow2 > 0 ) minpow2 = 0;
3254  }
3255  else {
3256  if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3257  if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3258  }
3259  else {
3260  MesPrint("Illegal term in expanded polyratfun.");
3261  goto onerror;
3262  }
3263  }
3264  rrm += *rrm;
3265  }
3266  AR.PolyFunPow += minpow2-1;
3267  }
3268  while ( m < arg2 ) {
3269  r = outarg;
3270  rrr = r++; mm = m + *m;
3271  *r++ = SYMBOL; *r++ = 4; *r++ = AR.PolyFunVar; *r++ = -1;
3272  m++; while ( m < mm ) *r++ = *m++;
3273  *rrr = r-rrr;
3274  Normalize(BHEAD rrr);
3275  StoreTerm(BHEAD rrr);
3276  }
3277  EndSort(BHEAD rr,1);
3278  r = rr; while ( *r ) r += *r;
3279  }
3280  if ( *rr == 0 ) {
3281  fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = CLEANFLAG;
3282  fun[1] = FUNHEAD+2;
3283  }
3284  else {
3285  rr = fun+FUNHEAD;
3286  *rr = r-rr;
3287  rr[1] = CLEANFLAG;
3288  if ( ToFast(rr,rr) ) {
3289  NEXTARG(rr);
3290  fun[1] = rr - fun;
3291  }
3292  else { fun[1] = r - fun; }
3293  }
3294  fun[2] = CLEANFLAG;
3295  goto done;
3296  }
3297  else if ( *arg2 == -SNUMBER ) {
3298  rr = r = outarg;
3299  if ( arg2[1] == 0 ) { error = 1; goto onerror; }
3300  if ( *arg1 == -SNUMBER ) { /* Things may not be normalized */
3301  if ( arg1[1] == 0 ) { *r++ = 0; }
3302  else {
3303  co1[0] = ABS(arg1[1]); co1[1] = 1;
3304  co2[0] = 1; co2[1] = ABS(arg2[1]);
3305  MulRat(BHEAD co1,1,co2,1,co,&nco);
3306  *r++ = 4; *r++ = (WORD)(co[0]); *r++ = (WORD)(co[1]);
3307  if ( ( arg1[1] < 0 && arg2[1] > 0 ) ||
3308  ( arg1[1] > 0 && arg2[1] < 0 ) ) *r++ = -3;
3309  else *r++ = 3;
3310  *r = 0;
3311  }
3312  }
3313  else if ( *arg1 == -SYMBOL ) {
3314  *r++ = 8; *r++ = SYMBOL; *r++ = 4;
3315  *r++ = arg1[1]; *r++ = 1;
3316  *r++ = 1; *r++ = ABS(arg2[1]);
3317  if ( arg2[1] < 0 ) *r++ = -3;
3318  else *r++ = 3;
3319  *r = 0;
3320  }
3321  else if ( *arg1 < 0 ) { error = 2; goto onerror; }
3322  else { /* Multi-term numerator. */
3323  m = arg1+ARGHEAD;
3324  NewSort(BHEAD0); /* Technically maybe not needed */
3325  if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3326  WORD minpow2 = MAXPOWER, *rrm;
3327  rrm = m;
3328  while ( rrm < arg2 ) {
3329  if ( *rrm == 4 ) {
3330  if ( minpow2 > 0 ) minpow2 = 0;
3331  }
3332  else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3333  if ( minpow2 > 0 ) minpow2 = 0;
3334  }
3335  else {
3336  if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3337  if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3338  }
3339  else {
3340  MesPrint("Illegal term in expanded polyratfun.");
3341  goto onerror;
3342  }
3343  }
3344  rrm += *rrm;
3345  }
3346  AR.PolyFunPow += minpow2;
3347  }
3348  while ( m < arg2 ) {
3349  r = rr;
3350  rrr = r++; mm = m + *m;
3351  *r++ = DENOMINATOR; *r++ = FUNHEAD + 2; *r++ = DIRTYFLAG;
3352  FILLFUN3(r);
3353  *r++ = arg2[0]; *r++ = arg2[1];
3354  m++; while ( m < mm ) *r++ = *m++;
3355  *rrr = r-rrr;
3356  if ( r < AT.WorkTop && r >= AT.WorkSpace )
3357  AT.WorkPointer = r;
3358  Normalize(BHEAD rrr);
3359  if ( ABS(rrr[*rrr-1]) == *rrr-1 ) {
3360  if ( AR.PolyFunPow >= 0 ) {
3361  StoreTerm(BHEAD rrr);
3362  }
3363  }
3364  else if ( rrr[1] == SYMBOL && rrr[2] == 4 &&
3365  rrr[3] == AR.PolyFunVar && rrr[4] <= AR.PolyFunPow ) {
3366  StoreTerm(BHEAD rrr);
3367  }
3368  }
3369  EndSort(BHEAD rr,1);
3370  }
3371  r = rr; while ( *r ) r += *r;
3372  i = r-rr;
3373  r = fun + FUNHEAD + ARGHEAD;
3374  NCOPY(r,rr,i);
3375  rr = fun + FUNHEAD;
3376  *rr = r - rr; rr[1] = CLEANFLAG;
3377  if ( ToFast(rr,rr) ) {
3378  NEXTARG(rr);
3379  fun[1] = rr - fun;
3380  }
3381  else { fun[1] = r - fun; }
3382  fun[2] = CLEANFLAG;
3383  goto done;
3384  }
3385  else { error = 0; goto onerror; }
3386  }
3387  else {
3388  r = arg2+ARGHEAD; /* The argument ends at tnext */
3389  first = 1;
3390  while ( r < tnext ) {
3391  rr = r + *r; rr -= ABS(rr[-1]);
3392  if ( r+1 == rr ) {
3393  if ( first ) { minpow = 0; first = 0; rmin = r; }
3394  else if ( minpow > 0 ) { minpow = 0; rmin = r; }
3395  }
3396  else if ( r[1] != SYMBOL || r[2] != 4 || r[3] != AR.PolyFunVar
3397  || r[4] > MAXPOWER ) { error = 0; goto onerror; }
3398  else if ( first ) { minpow = r[4]; first = 0; rmin = r; }
3399  else if ( r[4] < minpow ) { minpow = r[4]; rmin = r; }
3400  r += *r;
3401  }
3402 /*
3403  We have now:
3404  1: a numerator in arg1 which can contain several variables.
3405  2: a denominator in arg2 with at most only AR.PolyFunVar (ep).
3406  3: the minimum power in the denominator is minpow and the
3407  term with that minimum power is in rmin.
3408  Divide numerator and denominator by this minimum power.
3409  Determine the power range in the numerator.
3410  Call InvPoly.
3411  Multiply by the inverse in such a way that we never take more
3412  powers of ep than necessary.
3413 */
3414 /*
3415  One: put 1/rmin in AT.WorkPointer -> rmininv
3416 */
3417  AT.WorkPointer += AM.MaxTer/sizeof(WORD);
3418  if ( AT.WorkPointer + (AM.MaxTer/sizeof(WORD)) >= AT.WorkTop ) {
3419  error = 4; goto onerror;
3420  }
3421  rmininv = r = AT.WorkPointer;
3422  rr = rmin; i = *rmin; NCOPY(r,rr,i)
3423  if ( minpow != 0 ) { rmininv[4] = -rmininv[4]; }
3424  rsize = ABS(r[-1]);
3425  rcoef = r - rsize;
3426  rsize = (rsize-1)/2; rr = rcoef + rsize;
3427  for ( i = 0; i < rsize; i++ ) {
3428  rcopy = rcoef[i]; rcoef[i] = rr[i]; rr[i] = rcopy;
3429  }
3430  AT.WorkPointer = r;
3431  if ( *arg1 < 0 ) {
3432  ToGeneral(arg1,r,0);
3433  arg1 = r; r += *r; *r++ = 0; rcopy = 0;
3434  AT.WorkPointer = r;
3435  }
3436  else {
3437  r = arg1 + *arg1;
3438  rcopy = *r; *r++ = 0;
3439  }
3440 /*
3441  We can use MultiplyWithTerm.
3442 */
3443  AT.LeaveNegative = 1;
3444  numerator = MultiplyWithTerm(BHEAD arg1+ARGHEAD,rmininv,0);
3445  AT.LeaveNegative = 0;
3446  r[-1] = rcopy;
3447  r = numerator; while ( *r ) r += *r;
3448  AT.WorkPointer = r+1;
3449  rcopy = arg2[*arg2]; arg2[*arg2] = 0;
3450  denominator = MultiplyWithTerm(BHEAD arg2+ARGHEAD,rmininv,1);
3451  arg2[*arg2] = rcopy;
3452  r = denominator; while ( *r ) r += *r;
3453  AT.WorkPointer = r+1;
3454 /*
3455  Now find the minimum power of ep in the numerator.
3456 */
3457  r = numerator;
3458  first = 1;
3459  while ( *r ) {
3460  rr = r + *r; rr -= ABS(rr[-1]);
3461  if ( r+1 == rr ) {
3462  if ( first ) { minpow = 0; first = 0; }
3463  else if ( minpow > 0 ) { minpow = 0; }
3464  }
3465  else if ( r[1] != SYMBOL ) { error = 0; goto onerror; }
3466  else {
3467  for ( i = 3; i < r[2]; i += 2 ) {
3468  if ( r[i] == AR.PolyFunVar ) {
3469  if ( first ) { minpow = r[i+1]; first = 0; }
3470  else if ( r[i+1] < minpow ) minpow = r[i+1];
3471  break;
3472  }
3473  }
3474  if ( i >= r[2] ) {
3475  if ( first ) { minpow = 0; first = 0; }
3476  else if ( minpow > 0 ) minpow = 0;
3477  }
3478  }
3479  r += *r;
3480  }
3481 /*
3482  We can invert the denominator.
3483  Note that the return value is an offset in AT.pWorkSpace.
3484  Hence there is no need to free memory afterwards.
3485 */
3486  if ( AR.PolyFunExp == 3 ) {
3487  ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow-minpow,AR.PolyFunVar);
3488  }
3489  else {
3490  ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow,AR.PolyFunVar);
3491  }
3492 /*
3493  Now we start the multiplying
3494 */
3495  NewSort(BHEAD0);
3496  r = numerator;
3497  while ( *r ) {
3498 /*
3499  1: Find power of ep.
3500 */
3501  rnext = r + *r;
3502  rrr = rnext - ABS(rnext[-1]);
3503  rr = r+1;
3504  eppow = 0;
3505  if ( rr < rrr ) {
3506  j = rr[1] - 2; rr += 2;
3507  while ( j > 0 ) {
3508  if ( *rr == AR.PolyFunVar ) { eppow = rr[1]; break; }
3509  j -= 2; rr += 2;
3510  }
3511  }
3512 /*
3513  2: Multiply by the proper terms in ipoly
3514 */
3515  for ( i = 0; i <= AR.PolyFunPow-eppow+minpow; i++ ) {
3516  if ( AT.pWorkSpace[ipoly+i] == 0 ) continue;
3517 /*
3518  Copy the term, add i to the power of ep and multiply coef.
3519 */
3520  rc = r;
3521  rr = thecopy = AT.WorkPointer;
3522  while ( rc < rrr ) *rr++ = *rc++;
3523  if ( i != 0 ) {
3524  *rr++ = SYMBOL; *rr++ = 4; *rr++ = AR.PolyFunVar; *rr++ = i;
3525  }
3526  ncoef = REDLENG(rnext[-1]);
3527  MulRat(BHEAD (UWORD *)rrr,ncoef,
3528  (UWORD *)(AT.pWorkSpace[ipoly+i])+1,AT.pWorkSpace[ipoly+i][0]
3529  ,(UWORD *)rr,&newcoef);
3530  ncoef = ABS(newcoef); rr += 2*ncoef;
3531  newcoef = INCLENG(newcoef);
3532  *rr++ = newcoef;
3533  *thecopy = rr - thecopy;
3534  AT.WorkPointer = rr;
3535  Normalize(BHEAD thecopy);
3536  if ( *thecopy > 0 ) StoreTerm(BHEAD thecopy);
3537  AT.WorkPointer = thecopy;
3538  }
3539  r = rnext;
3540  }
3541 /*
3542  Now we have all.
3543 */
3544  rr = fun + FUNHEAD; r = rr + ARGHEAD;
3545  EndSort(BHEAD r,1);
3546  if ( *r == 0 ) {
3547  fun[1] = FUNHEAD+2; fun[2] = CLEANFLAG;
3548  fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0;
3549  }
3550  else {
3551  while ( *r ) r += *r;
3552  rr[0] = r-rr; rr[1] = CLEANFLAG;
3553  if ( ToFast(rr,rr) ) { NEXTARG(rr); fun[1] = rr-fun; }
3554  else { fun[1] = r-fun; }
3555  fun[2] = CLEANFLAG;
3556  }
3557  }
3558 done:
3559  if ( outarg ) TermFree(outarg-ARGHEAD,"ExpandRat");
3560  AR.PolyFunPow = OldPolyFunPow;
3561  AT.WorkPointer = ow;
3562  AN.PolyNormFlag = 1;
3563  return(0);
3564 onerror:
3565  if ( outarg ) TermFree(outarg-ARGHEAD,"ExpandRat");
3566  AR.PolyFunPow = OldPolyFunPow;
3567  AT.WorkPointer = ow;
3568  MLOCK(ErrorMessageLock);
3569  MesPrint(TheErrorMessage[error]);
3570  MUNLOCK(ErrorMessageLock);
3571  Terminate(-1);
3572  return(-1);
3573 }
3574 
3575 /*
3576  #] ExpandRat :
3577  #[ InvPoly :
3578 
3579  The input polynomial is represented as a sequence of terms in ascending
3580  power. The first coefficient is 1. If we call this 1-a and
3581  a = sum_(j,1,n,x^j*a(j)), and b = 1/(1-a) we can find the coefficients
3582  of b with the recursion
3583  b(0) = 1, b(n) = sum_(j,1,n,a(j)*b(n-j))
3584  The variable is the symbol sym and we need maxpow powers in the answer.
3585  The answer is an array of pointers to the coefficients of the various
3586  powers as rational numbers in the notation signedsize,numerator,denominator
3587  We put these powers in the workspace and the answer is in AT.pWorkSpace.
3588  Hence the return value is an offset in the pWorkSpace.
3589  A zero pointer indicates that this coefficient is zero.
3590 */
3591 
3592 int InvPoly(PHEAD WORD *inpoly, WORD maxpow, WORD sym)
3593 {
3594  int needed, inpointers, outpointers, maxinput = 0, i, j;
3595  WORD *t, *tt, *ttt, *w, *c, *cc, *ccc, lenc, lenc1, lenc2, rc, *c1, *c2;
3596 /*
3597  Step 0: allocate the space
3598 */
3599  needed = (maxpow+1)*2;
3600  WantAddPointers(needed);
3601  inpointers = AT.pWorkPointer;
3602  outpointers = AT.pWorkPointer+maxpow+1;
3603  for ( i = 0; i < needed; i++ ) AT.pWorkSpace[inpointers+i] = 0;
3604 /*
3605  Step 1: determine the coefficients in inpoly
3606  often there is a maximum power that is much smaller than maxpow.
3607  keeping track of this can speed up things.
3608 */
3609  t = inpoly;
3610  w = AT.WorkPointer;
3611  while ( *t ) {
3612  if ( *t == 4 ) {
3613  if ( t[1] != 1 || t[2] != 1 || t[3] != 3 ) goto onerror;
3614  AT.pWorkSpace[inpointers] = 0;
3615  }
3616  else if ( t[1] != SYMBOL || t[2] != 4 || t[3] != sym || t[4] < 0 ) goto onerror;
3617  else if ( t[4] > maxpow ) {} /* power outside useful range */
3618  else {
3619  if ( t[4] > maxinput ) maxinput = t[4];
3620  AT.pWorkSpace[inpointers+t[4]] = w;
3621  tt = t + *t; rc = -*--tt; /* we need - the coefficient! */
3622  rc = REDLENG(rc); *w++ = rc;
3623  ttt = t+5;
3624  while ( ttt < tt ) *w++ = *ttt++;
3625  }
3626  t += *t;
3627  }
3628 /*
3629  Step 2: compute the output. b(0) = 1.
3630  then the recursion starts.
3631 */
3632  AT.pWorkSpace[outpointers] = w;
3633  *w++ = 1; *w++ = 1; *w++ = 1;
3634  c = TermMalloc("InvPoly");
3635  c1 = TermMalloc("InvPoly");
3636  c2 = TermMalloc("InvPoly");
3637  for ( j = 1; j <= maxpow; j++ ) {
3638 /*
3639  Start at c = a(j)*b(0) = a(j)
3640 */
3641  if ( ( cc = AT.pWorkSpace[inpointers+j] ) != 0 ) {
3642  lenc = *cc++; /* reduced length */
3643  i = 2*ABS(lenc); ccc = c;
3644  NCOPY(ccc,cc,i);
3645  }
3646  else { lenc = 0; }
3647  for ( i = MiN(j-1,maxinput); i > 0; i-- ) {
3648 /*
3649  c -> c + a(i)*b(j-i)
3650 */
3651  if ( AT.pWorkSpace[inpointers+i] == 0
3652  || AT.pWorkSpace[outpointers+j-i] == 0 ) {
3653  }
3654  else {
3655  if ( MulRat(BHEAD (UWORD *)(AT.pWorkSpace[inpointers+i]+1),AT.pWorkSpace[inpointers+i][0],
3656  (UWORD *)(AT.pWorkSpace[outpointers+j-i]+1),AT.pWorkSpace[outpointers+j-i][0],
3657  (UWORD *)c1,&lenc1) ) goto calcerror;
3658  if ( lenc == 0 ) {
3659  cc = c; c = c1; c1 = cc;
3660  lenc = lenc1;
3661  }
3662  else {
3663  if ( AddRat(BHEAD (UWORD *)c,lenc,(UWORD *)c1,lenc1,(UWORD *)c2,&lenc2) )
3664  goto calcerror;
3665  cc = c; c = c2; c2 = cc;
3666  lenc = lenc2;
3667  }
3668  }
3669  }
3670 /*
3671  Copy c to the proper location
3672 */
3673  if ( lenc == 0 ) AT.pWorkSpace[outpointers+j] = 0;
3674  else {
3675  AT.pWorkSpace[outpointers+j] = w;
3676  *w++ = lenc;
3677  i = 2*ABS(lenc); ccc = c;
3678  NCOPY(w,ccc,i);
3679  }
3680  }
3681  AT.WorkPointer = w;
3682  TermFree(c2,"InvPoly");
3683  TermFree(c1,"InvPoly");
3684  TermFree(c ,"InvPoly");
3685 
3686  return(outpointers);
3687 onerror:
3688  MLOCK(ErrorMessageLock);
3689  MesPrint("Incorrect symbol field in InvPoly.");
3690  MUNLOCK(ErrorMessageLock);
3691  Terminate(-1);
3692  return(-1);
3693 calcerror:
3694  MLOCK(ErrorMessageLock);
3695  MesPrint("Called from InvPoly.");
3696  MUNLOCK(ErrorMessageLock);
3697  Terminate(-1);
3698  return(-1);
3699 }
3700 
3701 /*
3702  #] InvPoly :
3703 */
WORD CompareSymbols(WORD *, WORD *, WORD)
Definition: sort.c:2975
WORD * TakeSymbolContent(PHEAD WORD *, WORD *)
Definition: ratio.c:2434
Definition: structs.h:633
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
Definition: notation.c:510
int SymbolNormalize(WORD *)
Definition: normal.c:5000
Definition: structs.h:938
WORD InsertTerm(PHEAD WORD *, WORD, WORD, WORD *, WORD *, WORD)
Definition: proces.c:2550
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition: reken.c:1466
WORD StoreTerm(PHEAD WORD *)
Definition: sort.c:4332
WORD ** rhs
Definition: structs.h:943
VOID LowerSortLevel()
Definition: sort.c:4726
WORD * Buffer
Definition: structs.h:939
WORD NewSort(PHEAD0)
Definition: sort.c:591
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:3072
WORD * TakeContent(PHEAD WORD *, WORD *)
Definition: ratio.c:1376
LONG EndSort(PHEAD WORD *, int)
Definition: sort.c:681