FORM  4.2
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))<<1;
493  t += i++;
494  *t++ = (nc3 >= 0)?i:-i;
495  *term1 = WORDDIF(t,term1);
496 /*
497  Switch the new piece with the old tail, so that noncommuting
498  variables get into their proper spot.
499 */
500  i = old - insert;
501  t1 = t;
502  t2 = term1+insert;
503  NCOPY(t1,t2,i);
504  i = newer - old;
505  t1 = term1+insert;
506  t2 = term1+old;
507  NCOPY(t1,t2,i);
508  t2 = t;
509  i = old - insert;
510  NCOPY(t1,t2,i);
511  TermFree(coef,"Glue");
512  return(0);
513 }
514 
515 /*
516  #] Glue :
517  #[ DoSumF2 :
518 */
519 
520 WORD DoSumF2(PHEAD WORD *term, WORD *params, WORD replac, WORD level)
521 {
522  GETBIDENTITY
523  WORD *termout, *t, *from, *sub, *to, extractbuff = AT.TMbuff;
524  WORD isum, ival, iinc, insert, i;
525  CBUF *C;
526  ival = params[3];
527  iinc = params[5];
528  if ( ( iinc > 0 && params[4] >= ival )
529  || ( iinc < 0 && params[4] <= ival ) ) {
530  isum = (params[4] - ival)/iinc + 1;
531  }
532  else return(0);
533  termout = AT.WorkPointer;
534  AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
535  if ( AT.WorkPointer > AT.WorkTop ) {
536  MLOCK(ErrorMessageLock);
537  MesWork();
538  MUNLOCK(ErrorMessageLock);
539  return(-1);
540  }
541  t = term + 1;
542  while ( *t != SUBEXPRESSION || t[2] != replac || t[4] != extractbuff ) t += t[1];
543  insert = WORDDIF(t,term);
544 
545  from = term;
546  to = termout;
547  while ( from < t ) *to++ = *from++;
548  from += t[1];
549  sub = term + *term;
550  while ( from < sub ) *to++ = *from++;
551  *termout -= t[1];
552 
553  sub = t;
554  C = cbuf+t[4];
555  t += SUBEXPSIZE;
556  if ( params[2] < 0 ) {
557  while ( *t != INDTOIND || t[2] != -params[2] ) t += t[1];
558  *t = INDTOIND;
559  }
560  else {
561  while ( *t > SYMTOSUB || t[2] != params[2] ) t += t[1];
562  *t = SYMTONUM;
563  }
564  t[3] = ival;
565  for(;;) {
566  AT.WorkPointer = termout + *termout;
567  to = AT.WorkPointer;
568  if ( ( to + *termout ) > AT.WorkTop ) {
569  MLOCK(ErrorMessageLock);
570  MesWork();
571  MUNLOCK(ErrorMessageLock);
572  return(-1);
573  }
574  from = termout;
575  i = *termout;
576  NCOPY(to,from,i);
577  from = AT.WorkPointer;
578  AT.WorkPointer = to;
579  if ( Generator(BHEAD from,level) < 0 ) goto SumF2Call;
580  if ( --isum <= 0 ) break;
581  ival += iinc;
582  t[3] = ival;
583  if ( Glue(BHEAD termout,C->rhs[replac],sub,insert) < 0 ) goto SumF2Call;
584  }
585  AT.WorkPointer = termout;
586  return(0);
587 SumF2Call:
588  MLOCK(ErrorMessageLock);
589  MesCall("DoSumF2");
590  MUNLOCK(ErrorMessageLock);
591  SETERROR(-1)
592 }
593 
594 /*
595  #] DoSumF2 :
596  #] Sum :
597  #[ GCDfunction :
598  #[ GCDfunction :
599 */
600 
601 typedef struct {
602  WORD *buffer;
603  DOLLARS dollar;
604  LONG size;
605  int type;
606  int dummy;
607 } ARGBUFFER;
608 
609 int GCDfunction(PHEAD WORD *term,WORD level)
610 {
611  GETBIDENTITY
612  WORD *t, *tstop, *tf, *termout, *tin, *tout, *m, *mnext, *mstop, *mm;
613  int todo, i, ii, j, istart, sign = 1, action = 0;
614  WORD firstshort = 0, firstvalue = 0, gcdisone = 0, mlength, tlength, newlength;
615  WORD totargs = 0, numargs, 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  void *oldcompareroutine = AR.CompareRoutine;
1318  AR.CompareRoutine = (void *)&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 == 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 == 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 == 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 == 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 == 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  void *oldcompareroutine = AR.CompareRoutine;
2273  AR.SortType = SORTHIGHFIRST;
2274  AR.CompareRoutine = (void *)&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;
2794  int numargs = 0, type1, type2, actionflag1, actionflag2;
2795  WORD startebuf = cbuf[AT.ebufnum].numrhs;
2796  if ( par < 0 || par > 3 ) {
2797  MLOCK(ErrorMessageLock);
2798  MesPrint("Internal error. Illegal parameter %d in DIVfunction.",par);
2799  MUNLOCK(ErrorMessageLock);
2800  Terminate(-1);
2801  }
2802 /*
2803  Find the function
2804 */
2805  tend = term + *term; tstop = tend - ABS(tend[-1]);
2806  t = term+1;
2807  while ( t < tstop ) {
2808  if ( *t != divrem[par] ) { t += t[1]; continue; }
2809  r = t + FUNHEAD;
2810  tt = t + t[1]; numargs = 0;
2811  while ( r < tt ) {
2812  if ( numargs == 0 ) { arg1 = r; }
2813  if ( numargs == 1 ) { arg2 = r; }
2814  numargs++;
2815  NEXTARG(r);
2816  }
2817  if ( numargs == 2 ) break;
2818  t = tt;
2819  }
2820  if ( t >= tstop ) {
2821  MLOCK(ErrorMessageLock);
2822  MesPrint("Internal error. Indicated div_ or rem_ function not encountered.");
2823  MUNLOCK(ErrorMessageLock);
2824  Terminate(-1);
2825  }
2826 /*
2827  We have two arguments in arg1 and arg2.
2828 */
2829  if ( *arg1 == -SNUMBER && arg1[1] == 0 ) {
2830  if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2831 zerozero:;
2832  MLOCK(ErrorMessageLock);
2833  MesPrint("0/0 in either div_ or rem_ function.");
2834  MUNLOCK(ErrorMessageLock);
2835  Terminate(-1);
2836  }
2837  return(0);
2838  }
2839  if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2840 divzero:;
2841  MLOCK(ErrorMessageLock);
2842  MesPrint("Division by zero in either div_ or rem_ function.");
2843  MUNLOCK(ErrorMessageLock);
2844  Terminate(-1);
2845  }
2846  if ( ( arg1 = ConvertArgument(BHEAD arg1, &type1) ) == 0 ) goto CalledFrom;
2847  if ( ( arg2 = ConvertArgument(BHEAD arg2, &type2) ) == 0 ) goto CalledFrom;
2848  if ( *arg1 == 0 ) {
2849  if ( *arg2 == 0 ) {
2850  M_free(arg2,"DIVfunction");
2851  M_free(arg1,"DIVfunction");
2852  goto zerozero;
2853  }
2854  M_free(arg2,"DIVfunction");
2855  M_free(arg1,"DIVfunction");
2856  return(0);
2857  }
2858  if ( *arg2 == 0 ) {
2859  M_free(arg2,"DIVfunction");
2860  M_free(arg1,"DIVfunction");
2861  goto divzero;
2862  }
2863  if ( ( proper1 = PutExtraSymbols(BHEAD arg1,startebuf,&actionflag1) ) == 0 ) goto CalledFrom;
2864  if ( ( proper2 = PutExtraSymbols(BHEAD arg2,startebuf,&actionflag2) ) == 0 ) goto CalledFrom;
2865 /*
2866  if ( type2 == 0 ) M_free(arg2,"DIVfunction");
2867  else {
2868  DOLLARS d = ((DOLLARS)arg2)-1;
2869  if ( d->factors ) M_free(d->factors,"Dollar factors");
2870  M_free(d,"Copy of dollar variable");
2871  }
2872 */
2873  M_free(arg2,"DIVfunction");
2874 /*
2875  if ( type1 == 0 ) M_free(arg1,"DIVfunction");
2876  else {
2877  DOLLARS d = ((DOLLARS)arg1)-1;
2878  if ( d->factors ) M_free(d->factors,"Dollar factors");
2879  M_free(d,"Copy of dollar variable");
2880  }
2881 */
2882  M_free(arg1,"DIVfunction");
2883  if ( par == 0 ) proper3 = poly_div(BHEAD proper1, proper2,0);
2884  else if ( par == 1 ) proper3 = poly_rem(BHEAD proper1, proper2,0);
2885  else if ( par == 2 ) proper3 = poly_inverse(BHEAD proper1, proper2);
2886  else if ( par == 3 ) proper3 = poly_mul(BHEAD proper1, proper2);
2887  if ( proper3 == 0 ) goto CalledFrom;
2888  if ( actionflag1 || actionflag2 ) {
2889  if ( ( arg3 = TakeExtraSymbols(BHEAD proper3,startebuf) ) == 0 ) goto CalledFrom;
2890  M_free(proper3,"DIVfunction");
2891  }
2892  else {
2893  arg3 = proper3;
2894  }
2895  M_free(proper2,"DIVfunction");
2896  M_free(proper1,"DIVfunction");
2897  cbuf[AT.ebufnum].numrhs = startebuf;
2898  if ( *arg3 ) {
2899  termout = AT.WorkPointer;
2900  tlength = tend[-1];
2901  tlength = REDLENG(tlength);
2902  r3 = arg3;
2903  while ( *r3 ) {
2904  tt = term + 1; rr = termout + 1;
2905  while ( tt < t ) *rr++ = *tt++;
2906  r = r3 + 1;
2907  r3 = r3 + *r3;
2908  rstop = r3 - ABS(r3[-1]);
2909  while ( r < rstop ) *rr++ = *r++;
2910  tt += t[1];
2911  while ( tt < tstop ) *rr++ = *tt++;
2912  rlength = r3[-1];
2913  rlength = REDLENG(rlength);
2914  if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)rstop,rlength,
2915  (UWORD *)rr,&newlength) < 0 ) goto CalledFrom;
2916  rlength = INCLENG(newlength);
2917  rr += ABS(rlength);
2918  rr[-1] = rlength;
2919  *termout = rr - termout;
2920  AT.WorkPointer = rr;
2921  if ( Generator(BHEAD termout,level) ) goto CalledFrom;
2922  }
2923  AT.WorkPointer = termout;
2924  }
2925  M_free(arg3,"DIVfunction");
2926  return(0);
2927 CalledFrom:
2928  MLOCK(ErrorMessageLock);
2929  MesCall("DIVfunction");
2930  MUNLOCK(ErrorMessageLock);
2931  SETERROR(-1)
2932 }
2933 
2934 /*
2935  #] DIVfunction :
2936  #[ MULfunc :
2937 
2938  Multiplies two polynomials and puts the results in TermMalloc space.
2939 */
2940 
2941 WORD *MULfunc(PHEAD WORD *p1, WORD *p2)
2942 {
2943  WORD *prod,size1,size2,size3,*t,*tfill,*ps1,*ps2,sign1,sign2, error, *p3;
2944  UWORD *num1, *num2, *num3;
2945  int i;
2946  WORD oldsorttype = AR.SortType;
2947  void *oldcompareroutine = AR.CompareRoutine;
2948  AR.SortType = SORTHIGHFIRST;
2949  AR.CompareRoutine = (void *)&CompareSymbols;
2950  num3 = NumberMalloc("MULfunc");
2951  prod = TermMalloc("MULfunc");
2952  NewSort(BHEAD0);
2953  while ( *p1 ) {
2954  ps1 = p1+*p1; num1 = (UWORD *)(ps1 - ABS(ps1[-1])); size1 = ps1[-1];
2955  if ( size1 < 0 ) { sign1 = -1; size1 = -size1; }
2956  else sign1 = 1;
2957  size1 = (size1-1)/2;
2958  p3 = p2;
2959  while ( *p3 ) {
2960  ps2 = p3+*p3; num2 = (UWORD *)(ps2 - ABS(ps2[-1])); size2 = ps2[-1];
2961  if ( size2 < 0 ) { sign2 = -1; size2 = -size2; }
2962  else sign2 = 1;
2963  size2 = (size2-1)/2;
2964  if ( MulLong(num1,size1,num2,size2,num3,&size3) ) {
2965  error = 1;
2966 CalledFrom:
2967  MLOCK(ErrorMessageLock);
2968  MesPrint(" Error %d",error);
2969  MesCall("MulFunc");
2970  MUNLOCK(ErrorMessageLock);
2971  Terminate(-1);
2972  }
2973  tfill = prod+1;
2974  t = p1+1; while ( t < (WORD *)num1 ) *tfill++ = *t++;
2975  t = p3+1; while ( t < (WORD *)num2 ) *tfill++ = *t++;
2976  t = (WORD *)num3;
2977  for ( i = 0; i < size3; i++ ) *tfill++ = *t++;
2978  *tfill++ = 1;
2979  for ( i = 1; i < size3; i++ ) *tfill++ = 0;
2980  *tfill++ = (2*size3+1)*sign1*sign2;
2981  prod[0] = tfill - prod;
2982  if ( SymbolNormalize(prod) ) { error = 2; goto CalledFrom; }
2983  if ( StoreTerm(BHEAD prod) ) { error = 3; goto CalledFrom; }
2984  p3 += *p3;
2985  }
2986  p1 += *p1;
2987  }
2988  NumberFree(num3,"MULfunc");
2989  EndSort(BHEAD prod,1);
2990  AR.CompareRoutine = oldcompareroutine;
2991  AR.SortType = oldsorttype;
2992  return(prod);
2993 }
2994 
2995 /*
2996  #] MULfunc :
2997  #[ ConvertArgument :
2998 
2999  Converts an argument to a general notation in allocated space.
3000 */
3001 
3002 WORD *ConvertArgument(PHEAD WORD *arg, int *type)
3003 {
3004  WORD *output, *t, *r;
3005  int i, size;
3006  if ( *arg > 0 ) {
3007  output = (WORD *)Malloc1((*arg)*sizeof(WORD),"ConvertArgument");
3008  i = *arg - ARGHEAD; t = arg + ARGHEAD; r = output;
3009  NCOPY(r,t,i);
3010  *r = 0; *type = 0;
3011  return(output);
3012  }
3013  if ( *arg == -EXPRESSION ) {
3014  *type = 0;
3015  return(CreateExpression(BHEAD arg[1]));
3016  }
3017  if ( *arg == -DOLLAREXPRESSION ) {
3018  DOLLARS d;
3019  *type = 1;
3020  d = DolToTerms(BHEAD arg[1]);
3021 /*
3022  The problem is that DolToTerms creates a copy of the dollar variable.
3023  If we just return d->where we create a memory leak. Hence we have to
3024  copy the contents of d->where to a new buffer
3025 */
3026  output = (WORD *)Malloc1((d->size+1)*sizeof(WORD),"Copy of dollar content");
3027  WCOPY(output,d->where,d->size+1);
3028  if ( d->factors ) { M_free(d->factors,"Dollar factors"); d->factors = 0; }
3029  M_free(d,"Copy of dollar variable");
3030  return(output);
3031  }
3032 #if ( FUNHEAD > 4 )
3033  size = FUNHEAD+5;
3034 #else
3035  size = 9;
3036 #endif
3037  output = (WORD *)Malloc1(size*sizeof(WORD),"ConvertArgument");
3038  switch(*arg) {
3039  case -SYMBOL:
3040  output[0] = 8; output[1] = SYMBOL; output[2] = 4; output[3] = arg[1];
3041  output[4] = 1; output[5] = 1; output[6] = 1; output[7] = 3;
3042  output[8] = 0;
3043  break;
3044  case -INDEX:
3045  case -VECTOR:
3046  case -MINVECTOR:
3047  output[0] = 7; output[1] = INDEX; output[2] = 3; output[3] = arg[1];
3048  output[4] = 1; output[5] = 1;
3049  if ( *arg == -MINVECTOR ) output[6] = -3;
3050  else output[6] = 3;
3051  output[7] = 0;
3052  break;
3053  case -SNUMBER:
3054  output[0] = 4;
3055  if ( arg[1] < 0 ) {
3056  output[1] = -arg[1]; output[2] = 1; output[3] = -3;
3057  }
3058  else {
3059  output[1] = arg[1]; output[2] = 1; output[3] = 3;
3060  }
3061  output[4] = 0;
3062  break;
3063  default:
3064  output[0] = FUNHEAD+4;
3065  output[1] = -*arg;
3066  output[2] = FUNHEAD;
3067  for ( i = 3; i <= FUNHEAD; i++ ) output[i] = 0;
3068  output[FUNHEAD+1] = 1;
3069  output[FUNHEAD+2] = 1;
3070  output[FUNHEAD+3] = 3;
3071  output[FUNHEAD+4] = 0;
3072  break;
3073  }
3074  *type = 0;
3075  return(output);
3076 }
3077 
3078 /*
3079  #] ConvertArgument :
3080  #[ ExpandRat :
3081 
3082  Expands the denominator of a PolyRatFun in the variable PolyFunVar.
3083  The output is a polyratfun with a single argument.
3084  In the case that there is a polyratfun with more than one argument
3085  or the dirtyflag is on, the argument(s) is/are normalized.
3086  The output overwrites the input.
3087 */
3088 
3089 char *TheErrorMessage[] = {
3090  "PolyRatFun not of a type that FORM will expand: incorrect variable inside."
3091  ,"Division by zero in PolyRatFun encountered in ExpandRat."
3092  ,"Irregular code in PolyRatFun encountered in ExpandRat."
3093  ,"Called from ExpandRat."
3094  ,"WorkSpace overflow. Change parameter WorkSpace in setup file?"
3095  };
3096 
3097 int ExpandRat(PHEAD WORD *fun)
3098 {
3099  WORD *r, *rr, *rrr, *tt, *tnext, *arg1, *arg2, *rmin = 0, *rmininv;
3100  WORD *rcoef, rsize, rcopy, *ow = AT.WorkPointer;
3101  WORD *numerator, *denominator, *rnext;
3102  WORD *thecopy, *rc, ncoef, newcoef, *m, *mm, nco, *outarg = 0;
3103  UWORD co[2], co1[2], co2[2];
3104  WORD OldPolyFunPow = AR.PolyFunPow;
3105  int i, j, minpow = 0, eppow, first, error = 0, ipoly;
3106  if ( fun[1] == FUNHEAD ) { return(0); }
3107  tnext = fun + fun[1];
3108  if ( fun[1] == fun[FUNHEAD]+FUNHEAD ) { /* Single argument */
3109  if ( fun[2] == 0 ) { goto done; }
3110 /*
3111  We have to normalize the argument. This could make it shorter.
3112 */
3113 NormArg:;
3114  if ( outarg == 0 ) outarg = TermMalloc("ExpandRat")+ARGHEAD;
3115  AT.TrimPower = 1;
3116  NewSort(BHEAD0);
3117  r = fun+FUNHEAD+ARGHEAD;
3118  if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3119  WORD minpow2 = MAXPOWER, *rrm;
3120  rrm = r;
3121  while ( rrm < tnext ) {
3122  if ( *rrm == 4 ) {
3123  if ( minpow2 > 0 ) minpow2 = 0;
3124  }
3125  else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3126  if ( minpow2 > 0 ) minpow2 = 0;
3127  }
3128  else {
3129  if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3130  if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3131  }
3132  else {
3133  MesPrint("Illegal term in expanded polyratfun.");
3134  goto onerror;
3135  }
3136  }
3137  rrm += *rrm;
3138  }
3139  AR.PolyFunPow += minpow2;
3140  }
3141  while ( r < tnext ) {
3142  rr = r + *r;
3143  i = *r; rrr = outarg; NCOPY(rrr,r,i);
3144  Normalize(BHEAD outarg);
3145  if ( *outarg > 0 ) StoreTerm(BHEAD outarg);
3146  }
3147  r = fun+FUNHEAD+ARGHEAD;
3148  EndSort(BHEAD r,1);
3149  AT.TrimPower = 0;
3150  if ( *r == 0 ) {
3151  fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0;
3152  fun[1] = FUNHEAD+2;
3153  }
3154  else {
3155  rr = fun+FUNHEAD;
3156  if ( ToFast(rr,rr) ) {
3157  NEXTARG(rr); fun[1] = rr - fun;
3158  }
3159  else {
3160  while ( *r ) r += *r;
3161  *rr = r-rr; rr[1] = CLEANFLAG;
3162  fun[1] = r - fun;
3163  }
3164  }
3165  fun[2] = CLEANFLAG;
3166  goto done;
3167  }
3168 /*
3169  First test whether we have only AR.PolyFunVar in the denominator
3170 */
3171  tt = fun + FUNHEAD;
3172  arg1 = arg2 = 0;
3173  if ( tt < tnext ) {
3174  arg1 = tt; NEXTARG(tt);
3175  if ( tt < tnext ) {
3176  arg2 = tt; NEXTARG(tt);
3177  if ( tt != tnext ) { arg1 = arg2 = 0; } /* more than two arguments */
3178  }
3179  }
3180  if ( arg2 == 0 ) {
3181  if ( *arg1 < 0 ) { fun[2] = CLEANFLAG; goto done; }
3182  if ( fun[2] == CLEANFLAG ) goto done;
3183  goto NormArg; /* Note: should not come here */
3184  }
3185 /*
3186  Produce the output argument in outarg
3187 */
3188  if ( outarg == 0 ) outarg = TermMalloc("ExpandRat")+ARGHEAD;
3189 
3190  if ( *arg2 <= 0 ) {
3191 /*
3192  These cases are trivial.
3193  We try as much as possible to write the output directly into the
3194  function. We just have to be extremely careful not to overwrite
3195  relevant information before we are finished with it.
3196 */
3197  if ( *arg2 == -SYMBOL && arg2[1] == AR.PolyFunVar ) {
3198  rr = r = fun+FUNHEAD+ARGHEAD;
3199  if ( *arg1 < 0 ) {
3200  if ( *arg1 == -SYMBOL ) {
3201  if ( arg1[1] == AR.PolyFunVar ) {
3202  *r++ = 4; *r++ = 1; *r++ = 1; *r++ = 3; *r = 0;
3203  }
3204  else {
3205  *r++ = 10; *r++ = SYMBOL; *r++ = 6;
3206  *r++ = arg1[1]; *r++ = 1;
3207  *r++ = AR.PolyFunVar; *r++ = -1;
3208  *r++ = 1; *r++ = 1; *r++ = 3; *r = 0;
3209  Normalize(BHEAD rr);
3210  }
3211  }
3212  else if ( *arg1 == -SNUMBER ) {
3213  nco = arg1[1];
3214  if ( nco == 0 ) { *r++ = 0; }
3215  else {
3216  *r++ = 8; *r++ = SYMBOL; *r++ = 4;
3217  *r++ = AR.PolyFunVar; *r++ = -1;
3218  *r++ = ABS(nco); *r++ = 1;
3219  if ( nco < 0 ) *r++ = -3;
3220  else *r++ = 3;
3221  *r = 0;
3222  }
3223  }
3224  else { error = 2; goto onerror; } /* should not happen! */
3225  }
3226  else { /* Multi-term numerator. */
3227  m = arg1+ARGHEAD;
3228  NewSort(BHEAD0); /* Technically maybe not needed */
3229  if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3230  WORD minpow2 = MAXPOWER, *rrm;
3231  rrm = m;
3232  while ( rrm < arg2 ) {
3233  if ( *rrm == 4 ) {
3234  if ( minpow2 > 0 ) minpow2 = 0;
3235  }
3236  else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3237  if ( minpow2 > 0 ) minpow2 = 0;
3238  }
3239  else {
3240  if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3241  if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3242  }
3243  else {
3244  MesPrint("Illegal term in expanded polyratfun.");
3245  goto onerror;
3246  }
3247  }
3248  rrm += *rrm;
3249  }
3250  AR.PolyFunPow += minpow2-1;
3251  }
3252  while ( m < arg2 ) {
3253  r = outarg;
3254  rrr = r++; mm = m + *m;
3255  *r++ = SYMBOL; *r++ = 4; *r++ = AR.PolyFunVar; *r++ = -1;
3256  m++; while ( m < mm ) *r++ = *m++;
3257  *rrr = r-rrr;
3258  Normalize(BHEAD rrr);
3259  StoreTerm(BHEAD rrr);
3260  }
3261  EndSort(BHEAD rr,1);
3262  r = rr; while ( *r ) r += *r;
3263  }
3264  if ( *rr == 0 ) {
3265  fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = CLEANFLAG;
3266  fun[1] = FUNHEAD+2;
3267  }
3268  else {
3269  rr = fun+FUNHEAD;
3270  *rr = r-rr;
3271  rr[1] = CLEANFLAG;
3272  if ( ToFast(rr,rr) ) {
3273  NEXTARG(rr);
3274  fun[1] = rr - fun;
3275  }
3276  else { fun[1] = r - fun; }
3277  }
3278  fun[2] = CLEANFLAG;
3279  goto done;
3280  }
3281  else if ( *arg2 == -SNUMBER ) {
3282  rr = r = outarg;
3283  if ( arg2[1] == 0 ) { error = 1; goto onerror; }
3284  if ( *arg1 == -SNUMBER ) { /* Things may not be normalized */
3285  if ( arg1[1] == 0 ) { *r++ = 0; }
3286  else {
3287  co1[0] = ABS(arg1[1]); co1[1] = 1;
3288  co2[0] = 1; co2[1] = ABS(arg2[1]);
3289  MulRat(BHEAD co1,1,co2,1,co,&nco);
3290  *r++ = 4; *r++ = (WORD)(co[0]); *r++ = (WORD)(co[1]);
3291  if ( ( arg1[1] < 0 && arg2[1] > 0 ) ||
3292  ( arg1[1] > 0 && arg2[1] < 0 ) ) *r++ = -3;
3293  else *r++ = 3;
3294  *r = 0;
3295  }
3296  }
3297  else if ( *arg1 == -SYMBOL ) {
3298  *r++ = 8; *r++ = SYMBOL; *r++ = 4;
3299  *r++ = arg1[1]; *r++ = 1;
3300  *r++ = 1; *r++ = ABS(arg2[1]);
3301  if ( arg2[1] < 0 ) *r++ = -3;
3302  else *r++ = 3;
3303  *r = 0;
3304  }
3305  else if ( *arg1 < 0 ) { error = 2; goto onerror; }
3306  else { /* Multi-term numerator. */
3307  m = arg1+ARGHEAD;
3308  NewSort(BHEAD0); /* Technically maybe not needed */
3309  if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3310  WORD minpow2 = MAXPOWER, *rrm;
3311  rrm = m;
3312  while ( rrm < arg2 ) {
3313  if ( *rrm == 4 ) {
3314  if ( minpow2 > 0 ) minpow2 = 0;
3315  }
3316  else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3317  if ( minpow2 > 0 ) minpow2 = 0;
3318  }
3319  else {
3320  if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3321  if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3322  }
3323  else {
3324  MesPrint("Illegal term in expanded polyratfun.");
3325  goto onerror;
3326  }
3327  }
3328  rrm += *rrm;
3329  }
3330  AR.PolyFunPow += minpow2;
3331  }
3332  while ( m < arg2 ) {
3333  r = rr;
3334  rrr = r++; mm = m + *m;
3335  *r++ = DENOMINATOR; *r++ = FUNHEAD + 2; *r++ = DIRTYFLAG;
3336  FILLFUN3(r);
3337  *r++ = arg2[0]; *r++ = arg2[1];
3338  m++; while ( m < mm ) *r++ = *m++;
3339  *rrr = r-rrr;
3340  if ( r < AT.WorkTop && r >= AT.WorkSpace )
3341  AT.WorkPointer = r;
3342  Normalize(BHEAD rrr);
3343  if ( ABS(rrr[*rrr-1]) == *rrr-1 ) {
3344  if ( AR.PolyFunPow >= 0 ) {
3345  StoreTerm(BHEAD rrr);
3346  }
3347  }
3348  else if ( rrr[1] == SYMBOL && rrr[2] == 4 &&
3349  rrr[3] == AR.PolyFunVar && rrr[4] <= AR.PolyFunPow ) {
3350  StoreTerm(BHEAD rrr);
3351  }
3352  }
3353  EndSort(BHEAD rr,1);
3354  }
3355  r = rr; while ( *r ) r += *r;
3356  i = r-rr;
3357  r = fun + FUNHEAD + ARGHEAD;
3358  NCOPY(r,rr,i);
3359  rr = fun + FUNHEAD;
3360  *rr = r - rr; rr[1] = CLEANFLAG;
3361  if ( ToFast(rr,rr) ) {
3362  NEXTARG(rr);
3363  fun[1] = rr - fun;
3364  }
3365  else { fun[1] = r - fun; }
3366  fun[2] = CLEANFLAG;
3367  goto done;
3368  }
3369  else { error = 0; goto onerror; }
3370  }
3371  else {
3372  r = arg2+ARGHEAD; /* The argument ends at tnext */
3373  first = 1;
3374  while ( r < tnext ) {
3375  rr = r + *r; rr -= ABS(rr[-1]);
3376  if ( r+1 == rr ) {
3377  if ( first ) { minpow = 0; first = 0; rmin = r; }
3378  else if ( minpow > 0 ) { minpow = 0; rmin = r; }
3379  }
3380  else if ( r[1] != SYMBOL || r[2] != 4 || r[3] != AR.PolyFunVar
3381  || r[4] > MAXPOWER ) { error = 0; goto onerror; }
3382  else if ( first ) { minpow = r[4]; first = 0; rmin = r; }
3383  else if ( r[4] < minpow ) { minpow = r[4]; rmin = r; }
3384  r += *r;
3385  }
3386 /*
3387  We have now:
3388  1: a numerator in arg1 which can contain several variables.
3389  2: a denominator in arg2 with at most only AR.PolyFunVar (ep).
3390  3: the minimum power in the denominator is minpow and the
3391  term with that minimum power is in rmin.
3392  Divide numerator and denominator by this minimum power.
3393  Determine the power range in the numerator.
3394  Call InvPoly.
3395  Multiply by the inverse in such a way that we never take more
3396  powers of ep than necessary.
3397 */
3398 /*
3399  One: put 1/rmin in AT.WorkPointer -> rmininv
3400 */
3401  AT.WorkPointer += AM.MaxTer/sizeof(WORD);
3402  if ( AT.WorkPointer + (AM.MaxTer/sizeof(WORD)) >= AT.WorkTop ) {
3403  error = 4; goto onerror;
3404  }
3405  rmininv = r = AT.WorkPointer;
3406  rr = rmin; i = *rmin; NCOPY(r,rr,i)
3407  if ( minpow != 0 ) { rmininv[4] = -rmininv[4]; }
3408  rsize = ABS(r[-1]);
3409  rcoef = r - rsize;
3410  rsize = (rsize-1)/2; rr = rcoef + rsize;
3411  for ( i = 0; i < rsize; i++ ) {
3412  rcopy = rcoef[i]; rcoef[i] = rr[i]; rr[i] = rcopy;
3413  }
3414  AT.WorkPointer = r;
3415  if ( *arg1 < 0 ) {
3416  ToGeneral(arg1,r,0);
3417  arg1 = r; r += *r; *r++ = 0; rcopy = 0;
3418  AT.WorkPointer = r;
3419  }
3420  else {
3421  r = arg1 + *arg1;
3422  rcopy = *r; *r++ = 0;
3423  }
3424 /*
3425  We can use MultiplyWithTerm.
3426 */
3427  AT.LeaveNegative = 1;
3428  numerator = MultiplyWithTerm(BHEAD arg1+ARGHEAD,rmininv,0);
3429  AT.LeaveNegative = 0;
3430  r[-1] = rcopy;
3431  r = numerator; while ( *r ) r += *r;
3432  AT.WorkPointer = r+1;
3433  rcopy = arg2[*arg2]; arg2[*arg2] = 0;
3434  denominator = MultiplyWithTerm(BHEAD arg2+ARGHEAD,rmininv,1);
3435  arg2[*arg2] = rcopy;
3436  r = denominator; while ( *r ) r += *r;
3437  AT.WorkPointer = r+1;
3438 /*
3439  Now find the minimum power of ep in the numerator.
3440 */
3441  r = numerator;
3442  first = 1;
3443  while ( *r ) {
3444  rr = r + *r; rr -= ABS(rr[-1]);
3445  if ( r+1 == rr ) {
3446  if ( first ) { minpow = 0; first = 0; }
3447  else if ( minpow > 0 ) { minpow = 0; }
3448  }
3449  else if ( r[1] != SYMBOL ) { error = 0; goto onerror; }
3450  else {
3451  for ( i = 3; i < r[2]; i += 2 ) {
3452  if ( r[i] == AR.PolyFunVar ) {
3453  if ( first ) { minpow = r[i+1]; first = 0; }
3454  else if ( r[i+1] < minpow ) minpow = r[i+1];
3455  break;
3456  }
3457  }
3458  if ( i >= r[2] ) {
3459  if ( first ) { minpow = 0; first = 0; }
3460  else if ( minpow > 0 ) minpow = 0;
3461  }
3462  }
3463  r += *r;
3464  }
3465 /*
3466  We can invert the denominator.
3467  Note that the return value is an offset in AT.pWorkSpace.
3468  Hence there is no need to free memory afterwards.
3469 */
3470  if ( AR.PolyFunExp == 3 ) {
3471  ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow-minpow,AR.PolyFunVar);
3472  }
3473  else {
3474  ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow,AR.PolyFunVar);
3475  }
3476 /*
3477  Now we start the multiplying
3478 */
3479  NewSort(BHEAD0);
3480  r = numerator;
3481  while ( *r ) {
3482 /*
3483  1: Find power of ep.
3484 */
3485  rnext = r + *r;
3486  rrr = rnext - ABS(rnext[-1]);
3487  rr = r+1;
3488  eppow = 0;
3489  if ( rr < rrr ) {
3490  j = rr[1] - 2; rr += 2;
3491  while ( j > 0 ) {
3492  if ( *rr == AR.PolyFunVar ) { eppow = rr[1]; break; }
3493  j -= 2; rr += 2;
3494  }
3495  }
3496 /*
3497  2: Multiply by the proper terms in ipoly
3498 */
3499  for ( i = 0; i <= AR.PolyFunPow-eppow+minpow; i++ ) {
3500  if ( AT.pWorkSpace[ipoly+i] == 0 ) continue;
3501 /*
3502  Copy the term, add i to the power of ep and multiply coef.
3503 */
3504  rc = r;
3505  rr = thecopy = AT.WorkPointer;
3506  while ( rc < rrr ) *rr++ = *rc++;
3507  if ( i != 0 ) {
3508  *rr++ = SYMBOL; *rr++ = 4; *rr++ = AR.PolyFunVar; *rr++ = i;
3509  }
3510  ncoef = REDLENG(rnext[-1]);
3511  MulRat(BHEAD (UWORD *)rrr,ncoef,
3512  (UWORD *)(AT.pWorkSpace[ipoly+i])+1,AT.pWorkSpace[ipoly+i][0]
3513  ,(UWORD *)rr,&newcoef);
3514  ncoef = ABS(newcoef); rr += 2*ncoef;
3515  newcoef = INCLENG(newcoef);
3516  *rr++ = newcoef;
3517  *thecopy = rr - thecopy;
3518  AT.WorkPointer = rr;
3519  Normalize(BHEAD thecopy);
3520  if ( *thecopy > 0 ) StoreTerm(BHEAD thecopy);
3521  AT.WorkPointer = thecopy;
3522  }
3523  r = rnext;
3524  }
3525 /*
3526  Now we have all.
3527 */
3528  rr = fun + FUNHEAD; r = rr + ARGHEAD;
3529  EndSort(BHEAD r,1);
3530  if ( *r == 0 ) {
3531  fun[1] = FUNHEAD+2; fun[2] = CLEANFLAG;
3532  fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0;
3533  }
3534  else {
3535  while ( *r ) r += *r;
3536  rr[0] = r-rr; rr[1] = CLEANFLAG;
3537  if ( ToFast(rr,rr) ) { NEXTARG(rr); fun[1] = rr-fun; }
3538  else { fun[1] = r-fun; }
3539  fun[2] = CLEANFLAG;
3540  }
3541  }
3542 done:
3543  if ( outarg ) TermFree(outarg-ARGHEAD,"ExpandRat");
3544  AR.PolyFunPow = OldPolyFunPow;
3545  AT.WorkPointer = ow;
3546  AN.PolyNormFlag = 1;
3547  return(0);
3548 onerror:
3549  if ( outarg ) TermFree(outarg-ARGHEAD,"ExpandRat");
3550  AR.PolyFunPow = OldPolyFunPow;
3551  AT.WorkPointer = ow;
3552  MLOCK(ErrorMessageLock);
3553  MesPrint(TheErrorMessage[error]);
3554  MUNLOCK(ErrorMessageLock);
3555  Terminate(-1);
3556  return(-1);
3557 }
3558 
3559 /*
3560  #] ExpandRat :
3561  #[ InvPoly :
3562 
3563  The input polynomial is represented as a sequence of terms in ascending
3564  power. The first coefficient is 1. If we call this 1-a and
3565  a = sum_(j,1,n,x^j*a(j)), and b = 1/(1-a) we can find the coefficients
3566  of b with the recursion
3567  b(0) = 1, b(n) = sum_(j,1,n,a(j)*b(n-j))
3568  The variable is the symbol sym and we need maxpow powers in the answer.
3569  The answer is an array of pointers to the coefficients of the various
3570  powers as rational numbers in the notation signedsize,numerator,denominator
3571  We put these powers in the workspace and the answer is in AT.pWorkSpace.
3572  Hence the return value is an offset in the pWorkSpace.
3573  A zero pointer indicates that this coefficient is zero.
3574 */
3575 
3576 int InvPoly(PHEAD WORD *inpoly, WORD maxpow, WORD sym)
3577 {
3578  int needed, inpointers, outpointers, maxinput = 0, i, j;
3579  WORD *t, *tt, *ttt, *w, *c, *cc, *ccc, lenc, lenc1, lenc2, rc, *c1, *c2;
3580 /*
3581  Step 0: allocate the space
3582 */
3583  needed = (maxpow+1)*2;
3584  WantAddPointers(needed);
3585  inpointers = AT.pWorkPointer;
3586  outpointers = AT.pWorkPointer+maxpow+1;
3587  for ( i = 0; i < needed; i++ ) AT.pWorkSpace[inpointers+i] = 0;
3588 /*
3589  Step 1: determine the coefficients in inpoly
3590  often there is a maximum power that is much smaller than maxpow.
3591  keeping track of this can speed up things.
3592 */
3593  t = inpoly;
3594  w = AT.WorkPointer;
3595  while ( *t ) {
3596  if ( *t == 4 ) {
3597  if ( t[1] != 1 || t[2] != 1 || t[3] != 3 ) goto onerror;
3598  AT.pWorkSpace[inpointers] = 0;
3599  }
3600  else if ( t[1] != SYMBOL || t[2] != 4 || t[3] != sym || t[4] < 0 ) goto onerror;
3601  else if ( t[4] > maxpow ) {} /* power outside useful range */
3602  else {
3603  if ( t[4] > maxinput ) maxinput = t[4];
3604  AT.pWorkSpace[inpointers+t[4]] = w;
3605  tt = t + *t; rc = -*--tt; /* we need - the coefficient! */
3606  rc = REDLENG(rc); *w++ = rc;
3607  ttt = t+5;
3608  while ( ttt < tt ) *w++ = *ttt++;
3609  }
3610  t += *t;
3611  }
3612 /*
3613  Step 2: compute the output. b(0) = 1.
3614  then the recursion starts.
3615 */
3616  AT.pWorkSpace[outpointers] = w;
3617  *w++ = 1; *w++ = 1; *w++ = 1;
3618  c = TermMalloc("InvPoly");
3619  c1 = TermMalloc("InvPoly");
3620  c2 = TermMalloc("InvPoly");
3621  for ( j = 1; j <= maxpow; j++ ) {
3622 /*
3623  Start at c = a(j)*b(0) = a(j)
3624 */
3625  if ( ( cc = AT.pWorkSpace[inpointers+j] ) != 0 ) {
3626  lenc = *cc++; /* reduced length */
3627  i = 2*ABS(lenc); ccc = c;
3628  NCOPY(ccc,cc,i);
3629  }
3630  else { lenc = 0; }
3631  for ( i = MiN(j-1,maxinput); i > 0; i-- ) {
3632 /*
3633  c -> c + a(i)*b(j-i)
3634 */
3635  if ( AT.pWorkSpace[inpointers+i] == 0
3636  || AT.pWorkSpace[outpointers+j-i] == 0 ) {
3637  }
3638  else {
3639  if ( MulRat(BHEAD (UWORD *)(AT.pWorkSpace[inpointers+i]+1),AT.pWorkSpace[inpointers+i][0],
3640  (UWORD *)(AT.pWorkSpace[outpointers+j-i]+1),AT.pWorkSpace[outpointers+j-i][0],
3641  (UWORD *)c1,&lenc1) ) goto calcerror;
3642  if ( lenc == 0 ) {
3643  cc = c; c = c1; c1 = cc;
3644  lenc = lenc1;
3645  }
3646  else {
3647  if ( AddRat(BHEAD (UWORD *)c,lenc,(UWORD *)c1,lenc1,(UWORD *)c2,&lenc2) )
3648  goto calcerror;
3649  cc = c; c = c2; c2 = cc;
3650  lenc = lenc2;
3651  }
3652  }
3653  }
3654 /*
3655  Copy c to the proper location
3656 */
3657  if ( lenc == 0 ) AT.pWorkSpace[outpointers+j] = 0;
3658  else {
3659  AT.pWorkSpace[outpointers+j] = w;
3660  *w++ = lenc;
3661  i = 2*ABS(lenc); ccc = c;
3662  NCOPY(w,ccc,i);
3663  }
3664  }
3665  AT.WorkPointer = w;
3666  TermFree(c2,"InvPoly");
3667  TermFree(c1,"InvPoly");
3668  TermFree(c ,"InvPoly");
3669 
3670  return(outpointers);
3671 onerror:
3672  MLOCK(ErrorMessageLock);
3673  MesPrint("Incorrect symbol field in InvPoly.");
3674  MUNLOCK(ErrorMessageLock);
3675  Terminate(-1);
3676  return(-1);
3677 calcerror:
3678  MLOCK(ErrorMessageLock);
3679  MesPrint("Called from InvPoly.");
3680  MUNLOCK(ErrorMessageLock);
3681  Terminate(-1);
3682  return(-1);
3683 }
3684 
3685 /*
3686  #] InvPoly :
3687 */
WORD * TakeSymbolContent(PHEAD WORD *in, WORD *term)
Definition: ratio.c:2434
WORD * TakeContent(PHEAD WORD *in, WORD *term)
Definition: ratio.c:1376
Definition: structs.h:620
#define PHEAD
Definition: ftypes.h:56
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
Definition: notation.c:510
int SymbolNormalize(WORD *)
Definition: normal.c:4979
Definition: structs.h:921
WORD InsertTerm(PHEAD WORD *, WORD, WORD, WORD *, WORD *, WORD)
Definition: proces.c:2514
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition: reken.c:1466
WORD StoreTerm(PHEAD WORD *)
Definition: sort.c:4246
WORD ** rhs
Definition: structs.h:926
VOID LowerSortLevel()
Definition: sort.c:4610
WORD * Buffer
Definition: structs.h:922
WORD NewSort(PHEAD0)
Definition: sort.c:589
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:3034
int CompareSymbols(PHEAD WORD *, WORD *, WORD)
Definition: sort.c:2945
LONG EndSort(PHEAD WORD *, int)
Definition: sort.c:675