FORM 4.3
ratio.c
Go to the documentation of this file.
1
8/* #[ License : */
9/*
10 * Copyright (C) 1984-2022 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
62WORD 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
93We 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
170WORD 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;
212PosNeg:
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
305RatioCall:
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
320WORD 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
395WORD 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);
443SumF1Call:
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
461WORD Glue(PHEAD WORD *term1, WORD *term2, WORD *sub, WORD insert)
462{
463 GETBIDENTITY
464 UWORD *coef;
465 WORD ncoef, *t, *t1, *t2, i, nc2, nc3, old, newer;
466 coef = (UWORD *)(TermMalloc("Glue"));
467 t = term1;
468 t += *t;
469 i = t[-1];
470 t -= ABS(i);
471 old = WORDDIF(t,term1);
472 ncoef = REDLENG(i);
473 if ( i < 0 ) i = -i;
474 i--;
475 t1 = t;
476 t2 = (WORD *)coef;
477 while ( --i >= 0 ) *t2++ = *t1++;
478 i = *--t;
479 nc2 = WildFill(BHEAD t,term2,sub);
480 *t = i;
481 t += nc2;
482 nc2 = t[-1];
483 t -= ABS(nc2);
484 newer = WORDDIF(t,term1);
485 if ( MulRat(BHEAD (UWORD *)t,REDLENG(nc2),coef,ncoef,(UWORD *)t,&nc3) ) {
486 MLOCK(ErrorMessageLock);
487 MesCall("Glue");
488 MUNLOCK(ErrorMessageLock);
489 TermFree(coef,"Glue");
490 SETERROR(-1)
491 }
492 i = (ABS(nc3))*2;
493 t += i++;
494 *t++ = (nc3 >= 0)?i:-i;
495 *term1 = WORDDIF(t,term1);
496/*
497 Switch the new piece with the old tail, so that noncommuting
498 variables get into their proper spot.
499*/
500 i = old - insert;
501 t1 = t;
502 t2 = term1+insert;
503 NCOPY(t1,t2,i);
504 i = newer - old;
505 t1 = term1+insert;
506 t2 = term1+old;
507 NCOPY(t1,t2,i);
508 t2 = t;
509 i = old - insert;
510 NCOPY(t1,t2,i);
511 TermFree(coef,"Glue");
512 return(0);
513}
514
515/*
516 #] Glue :
517 #[ DoSumF2 :
518*/
519
520WORD 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);
587SumF2Call:
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
601typedef struct {
602 WORD *buffer;
603 DOLLARS dollar;
604 LONG size;
605 int type;
606 int dummy;
607} ARGBUFFER;
608
609int 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;
733gcdone:
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 */
749doshort:
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;
852oneterm:;
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;
1068multiterms:;
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*/
1093cleanup:;
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
1113CalledFrom:
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
1141WORD *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);
1231CalledFrom:
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
1244WORD *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 ) {
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);
1261CalledFrom:
1262 MLOCK(ErrorMessageLock);
1263 MesCall("PutExtraSymbols");
1264 MUNLOCK(ErrorMessageLock);
1265 return(0);
1266}
1267
1268/*
1269 #] PutExtraSymbols :
1270 #[ TakeExtraSymbols :
1271*/
1272
1273WORD *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 ) {
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) ) {
1293 goto CalledFrom;
1294 }
1295 }
1296 AT.WorkPointer = oldworkpointer;
1297 if ( EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 ) goto CalledFrom;
1298 return(termout);
1299
1300CalledFrom:
1301 MLOCK(ErrorMessageLock);
1302 MesCall("TakeExtraSymbols");
1303 MUNLOCK(ErrorMessageLock);
1304 return(0);
1305}
1306
1307/*
1308 #] TakeExtraSymbols :
1309 #[ MultiplyWithTerm :
1310*/
1311
1312WORD *MultiplyWithTerm(PHEAD WORD *in, WORD *term, WORD par)
1313{
1314 WORD *termout, *t, *tt, *tstop, *ttstop;
1315 WORD length, length1, length2;
1316 WORD oldsorttype = AR.SortType;
1317 COMPARE oldcompareroutine = AR.CompareRoutine;
1318 AR.CompareRoutine = &CompareSymbols;
1319
1320 if ( par == 0 || par == 2 ) AR.SortType = SORTHIGHFIRST;
1321 else AR.SortType = SORTLOWFIRST;
1322 termout = AT.WorkPointer;
1323 NewSort(BHEAD0);
1324 while ( *in ) {
1325 tt = termout + 1;
1326 tstop = in + *in; tstop -= ABS(tstop[-1]); t = in + 1;
1327 while ( t < tstop ) *tt++ = *t++;
1328 ttstop = term + *term; ttstop -= ABS(ttstop[-1]); t = term + 1;
1329 while ( t < ttstop ) *tt++ = *t++;
1330 length1 = REDLENG(in[*in-1]); length2 = REDLENG(term[*term-1]);
1331 if ( MulRat(BHEAD (UWORD *)tstop,length1,
1332 (UWORD *)ttstop,length2,(UWORD *)tt,&length) ) goto CalledFrom;
1333 length = INCLENG(length);
1334 tt += ABS(length); tt[-1] = length;
1335 *termout = tt - termout;
1336 SymbolNormalize(termout);
1337 StoreTerm(BHEAD termout);
1338 in += *in;
1339 }
1340 if ( par == 2 ) {
1341/* if ( AN.tryterm == 0 ) AN.tryterm = 1; */
1342 AN.tryterm = 0; /* For now */
1343 if ( EndSort(BHEAD (WORD *)((VOID *)(&termout)),2) < 0 ) goto CalledFrom;
1344 }
1345 else {
1346 if ( EndSort(BHEAD termout,1) < 0 ) goto CalledFrom;
1347 }
1348
1349 AR.CompareRoutine = oldcompareroutine;
1350
1351 AR.SortType = oldsorttype;
1352 return(termout);
1353
1354CalledFrom:
1355 MLOCK(ErrorMessageLock);
1356 MesCall("MultiplyWithTerm");
1357 MUNLOCK(ErrorMessageLock);
1358 return(0);
1359}
1360
1361/*
1362 #] MultiplyWithTerm :
1363 #[ TakeContent :
1364*/
1375
1376WORD *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 }
1446noindices:
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 }
1516novectors:;
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;
1565nextr1:;
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 }
1600nofunctions:
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);
1789CalledFrom:
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
1808int 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
1927int 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
2020WORD *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) ) {
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);
2059CalledFrom:
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
2076int GCDterms(PHEAD WORD *term1, WORD *term2, WORD *termout)
2077{
2078 GETBIDENTITY
2079 WORD *t1, *t1stop, *t1next, *t2, *t2stop, *t2next, *tout, *tt1, *tt2;
2080 int count1, count2, i, ii, x1, sign;
2081 WORD length1, length2;
2082 t1 = term1 + *term1; t1stop = t1 - ABS(t1[-1]); t1 = term1+1;
2083 t2 = term2 + *term2; t2stop = t2 - ABS(t2[-1]); t2 = term2+1;
2084 tout = termout+1;
2085 while ( t1 < t1stop ) {
2086 t1next = t1 + t1[1];
2087 t2 = term2+1;
2088 if ( *t1 == SYMBOL ) {
2089 while ( t2 < t2stop && *t2 != SYMBOL ) t2 += t2[1];
2090 if ( t2 < t2stop && *t2 == SYMBOL ) {
2091 t2next = t2+t2[1];
2092 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2093 while ( tt1 < t1next && tt2 < t2next ) {
2094 if ( *tt1 < *tt2 ) tt1 += 2;
2095 else if ( *tt1 > *tt2 ) tt2 += 2;
2096 else if ( ( tt1[1] > 0 && tt2[1] < 0 ) ||
2097 ( tt2[1] > 0 && tt1[1] < 0 ) ) {
2098 tt1 += 2; tt2 += 2;
2099 }
2100 else {
2101 x1 = tt1[1];
2102 if ( tt1[1] < 0 ) { if ( tt2[1] > x1 ) x1 = tt2[1]; }
2103 else { if ( tt2[1] < x1 ) x1 = tt2[1]; }
2104 tout[count1+2] = *tt1;
2105 tout[count1+3] = x1;
2106 tt1 += 2; tt2 += 2;
2107 count1 += 2;
2108 }
2109 }
2110 if ( count1 > 0 ) {
2111 *tout = SYMBOL; tout[1] = count1+2; tout += tout[1];
2112 }
2113 }
2114 }
2115 else if ( *t1 == DOTPRODUCT ) {
2116 while ( t2 < t2stop && *t2 != DOTPRODUCT ) t2 += t2[1];
2117 if ( t2 < t2stop && *t2 == DOTPRODUCT ) {
2118 t2next = t2+t2[1];
2119 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2120 while ( tt1 < t1next && tt2 < t2next ) {
2121 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 3;
2122 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 3;
2123 else if ( ( tt1[2] > 0 && tt2[2] < 0 ) ||
2124 ( tt2[2] > 0 && tt1[2] < 0 ) ) {
2125 tt1 += 3; tt2 += 3;
2126 }
2127 else {
2128 x1 = tt1[2];
2129 if ( tt1[2] < 0 ) { if ( tt2[2] > x1 ) x1 = tt2[2]; }
2130 else { if ( tt2[2] < x1 ) x1 = tt2[2]; }
2131 tout[count1+2] = *tt1;
2132 tout[count1+3] = tt1[1];
2133 tout[count1+4] = x1;
2134 tt1 += 3; tt2 += 3;
2135 count1 += 3;
2136 }
2137 }
2138 if ( count1 > 0 ) {
2139 *tout = DOTPRODUCT; tout[1] = count1+2; tout += tout[1];
2140 }
2141 }
2142 }
2143 else if ( *t1 == VECTOR ) {
2144 while ( t2 < t2stop && *t2 != VECTOR ) t2 += t2[1];
2145 if ( t2 < t2stop && *t2 == VECTOR ) {
2146 t2next = t2+t2[1];
2147 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2148 while ( tt1 < t1next && tt2 < t2next ) {
2149 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2150 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2151 else {
2152 tout[count1+2] = *tt1;
2153 tout[count1+3] = tt1[1];
2154 tt1 += 2; tt2 += 2;
2155 count1 += 2;
2156 }
2157 }
2158 if ( count1 > 0 ) {
2159 *tout = VECTOR; tout[1] = count1+2; tout += tout[1];
2160 }
2161 }
2162 }
2163 else if ( *t1 == INDEX ) {
2164 while ( t2 < t2stop && *t2 != INDEX ) t2 += t2[1];
2165 if ( t2 < t2stop && *t2 == INDEX ) {
2166 t2next = t2+t2[1];
2167 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2168 while ( tt1 < t1next && tt2 < t2next ) {
2169 if ( *tt1 < *tt2 ) tt1 += 1;
2170 else if ( *tt1 > *tt2 ) tt2 += 1;
2171 else {
2172 tout[count1+2] = *tt1;
2173 tt1 += 1; tt2 += 1;
2174 count1 += 1;
2175 }
2176 }
2177 if ( count1 > 0 ) {
2178 *tout = INDEX; tout[1] = count1+2; tout += tout[1];
2179 }
2180 }
2181 }
2182 else if ( *t1 == DELTA ) {
2183 while ( t2 < t2stop && *t2 != DELTA ) t2 += t2[1];
2184 if ( t2 < t2stop && *t2 == DELTA ) {
2185 t2next = t2+t2[1];
2186 tt1 = t1+2; tt2 = t2+2; count1 = 0;
2187 while ( tt1 < t1next && tt2 < t2next ) {
2188 if ( *tt1 < *tt2 || ( *tt1 == *tt2 && tt1[1] < tt2[1] ) ) tt1 += 2;
2189 else if ( *tt1 > *tt2 || ( *tt1 == *tt2 && tt1[1] > tt2[1] ) ) tt2 += 2;
2190 else {
2191 tout[count1+2] = *tt1;
2192 tout[count1+3] = tt1[1];
2193 tt1 += 2; tt2 += 2;
2194 count1 += 2;
2195 }
2196 }
2197 if ( count1 > 0 ) {
2198 *tout = DELTA; tout[1] = count1+2; tout += tout[1];
2199 }
2200 }
2201 }
2202 else if ( *t1 >= FUNCTION ) { /* noncommuting functions? Forbidden! */
2203/*
2204 Count how many times this function occurs.
2205 Then count how many times it is in term2.
2206*/
2207 count1 = 1;
2208 while ( t1next < t1stop && *t1 == *t1next && t1[1] == t1next[1] ) {
2209 for ( i = 2; i < t1[1]; i++ ) {
2210 if ( t1[i] != t1next[i] ) break;
2211 }
2212 if ( i < t1[1] ) break;
2213 count1++;
2214 t1next += t1next[1];
2215 }
2216 count2 = 0;
2217 while ( t2 < t2stop ) {
2218 if ( *t2 == *t1 && t2[1] == t1[1] ) {
2219 for ( i = 2; i < t1[1]; i++ ) {
2220 if ( t2[i] != t1[i] ) break;
2221 }
2222 if ( i >= t1[1] ) count2++;
2223 }
2224 t2 += t2[1];
2225 }
2226 if ( count1 < count2 ) count2 = count1; /* number of common occurrences */
2227 if ( count2 > 0 ) {
2228 if ( tout == t1 ) {
2229 while ( count2 > 0 ) { tout += tout[1]; count2--; }
2230 }
2231 else {
2232 i = t1[1]*count2;
2233 NCOPY(tout,t1,i);
2234 }
2235 }
2236 }
2237 t1 = t1next;
2238 }
2239/*
2240 Now the coefficients. They are in t1stop and t2stop. Should go to tout.
2241*/
2242 sign = 1;
2243 length1 = term1[*term1-1]; ii = i = ABS(length1); t1 = t1stop;
2244 if ( t1 != tout ) { NCOPY(tout,t1,i); tout -= ii; }
2245 length2 = term2[*term2-1];
2246 if ( length1 < 0 && length2 < 0 ) sign = -1;
2247 if ( AccumGCD(BHEAD (UWORD *)tout,&length1,(UWORD *)t2stop,length2) ) {
2248 MLOCK(ErrorMessageLock);
2249 MesCall("GCDterms");
2250 MUNLOCK(ErrorMessageLock);
2251 SETERROR(-1)
2252 }
2253 if ( sign < 0 && length1 > 0 ) length1 = -length1;
2254 tout += ABS(length1); tout[-1] = length1;
2255 *termout = tout - termout; *tout = 0;
2256 return(0);
2257}
2258
2259/*
2260 #] GCDterms :
2261 #[ ReadPolyRatFun :
2262*/
2263
2264int ReadPolyRatFun(PHEAD WORD *term)
2265{
2266 WORD *oldworkpointer = AT.WorkPointer;
2267 int flag, i;
2268 WORD *t, *fun, *nextt, *num, *den, *t1, *t2, size, numsize, densize;
2269 WORD *term1, *term2, *confree1, *confree2, *gcd, *num1, *den1, move, *newnum, *newden;
2270 WORD *tstop, *m1, *m2;
2271 WORD oldsorttype = AR.SortType;
2272 COMPARE oldcompareroutine = AR.CompareRoutine;
2273 AR.SortType = SORTHIGHFIRST;
2274 AR.CompareRoutine = &CompareSymbols;
2275
2276 tstop = term + *term; tstop -= ABS(tstop[-1]);
2277 if ( term + *term == AT.WorkPointer ) flag = 1;
2278 else flag = 0;
2279 t = term+1;
2280 while ( t < tstop ) {
2281 if ( *t != AR.PolyFun ) { t += t[1]; continue; }
2282 if ( ( t[2] & MUSTCLEANPRF ) == 0 ) { t += t[1]; continue; }
2283 fun = t;
2284 nextt = t + t[1];
2285 if ( fun[1] > FUNHEAD && fun[FUNHEAD] == -SNUMBER && fun[FUNHEAD+1] == 0 )
2286 { *term = 0; break; }
2287 if ( FromPolyRatFun(BHEAD fun, &num, &den) > 0 ) { t = nextt; continue; }
2288 if ( *num == ARGHEAD ) { *term = 0; break; }
2289/*
2290 Now we have num and den. Both are in general argument notation,
2291 but can also be used as expressions as in num+ARGHEAD, den+ARGHEAD.
2292 We need the gcd. For this we have to take out the contents
2293 because PreGCD does not like contents.
2294*/
2295 term1 = TermMalloc("ReadPolyRatFun");
2296 term2 = TermMalloc("ReadPolyRatFun");
2297 confree1 = TakeSymbolContent(BHEAD num+ARGHEAD,term1);
2298 confree2 = TakeSymbolContent(BHEAD den+ARGHEAD,term2);
2299 GCDclean(BHEAD term1,term2);
2300/* gcd = PreGCD(BHEAD confree1,confree2,1); */
2301 gcd = poly_gcd(BHEAD confree1,confree2,1);
2302 newnum = PolyDiv(BHEAD confree1,gcd,"ReadPolyRatFun");
2303 newden = PolyDiv(BHEAD confree2,gcd,"ReadPolyRatFun");
2304 TermFree(confree2,"ReadPolyRatFun");
2305 TermFree(confree1,"ReadPolyRatFun");
2306 num1 = MULfunc(BHEAD term1,newnum);
2307 den1 = MULfunc(BHEAD term2,newden);
2308 TermFree(newnum,"ReadPolyRatFun");
2309 TermFree(newden,"ReadPolyRatFun");
2310/* M_free(gcd,"poly_gcd"); */
2311 TermFree(gcd,"poly_gcd");
2312 TermFree(term1,"ReadPolyRatFun");
2313 TermFree(term2,"ReadPolyRatFun");
2314/*
2315 Now we can put the function back together.
2316 Notice that we cannot use ToFast, because there is no reservation
2317 for the header of the argument. Fortunately there are only two
2318 types of fast arguments.
2319*/
2320 if ( num1[0] == 4 && num1[4] == 0 && num1[2] == 1 && num1[1] > 0 ) {
2321 numsize = 2; num1[0] = -SNUMBER;
2322 if ( num1[3] < 0 ) num1[1] = -num1[1];
2323 }
2324 else if ( num1[0] == 8 && num1[8] == 0 && num1[7] == 3 && num1[6] == 1
2325 && num1[5] == 1 && num1[1] == SYMBOL && num1[4] == 1 ) {
2326 numsize = 2; num1[0] = -SYMBOL; num1[1] = num1[3];
2327 }
2328 else { m1 = num1; while ( *m1 ) m1 += *m1; numsize = (m1-num1)+ARGHEAD; }
2329 if ( den1[0] == 4 && den1[4] == 0 && den1[2] == 1 && den1[1] > 0 ) {
2330 densize = 2; den1[0] = -SNUMBER;
2331 if ( den1[3] < 0 ) den1[1] = -den1[1];
2332 }
2333 else if ( den1[0] == 8 && den1[8] == 0 && den1[7] == 3 && den1[6] == 1
2334 && den1[5] == 1 && den1[1] == SYMBOL && den1[4] == 1 ) {
2335 densize = 2; den1[0] = -SYMBOL; den1[1] = den1[3];
2336 }
2337 else { m2 = den1; while ( *m2 ) m2 += *m2; densize = (m2-den1)+ARGHEAD; }
2338 size = FUNHEAD+numsize+densize;
2339
2340 if ( size > fun[1] ) {
2341 move = size - fun[1];
2342 t1 = term+*term; t2 = t1+move;
2343 while ( t1 > nextt ) *--t2 = *--t1;
2344 tstop += move; nextt += move;
2345 *term += move;
2346 }
2347 else if ( size < fun[1] ) {
2348 move = fun[1]-size;
2349 t2 = fun+size; t1 = nextt;
2350 tstop -= move; nextt -= move;
2351 t = term+*term;
2352 while ( t1 < t ) *t2++ = *t1++;
2353 *term -= move;
2354 }
2355 else { /* no need to move anything */ }
2356 fun[1] = size; fun[2] = 0;
2357 t2 = fun+FUNHEAD; t1 = num1;
2358 if ( *num1 < 0 ) { *t2++ = num1[0]; *t2++ = num1[1]; }
2359 else { *t2++ = numsize; *t2++ = 0; FILLARG(t2);
2360 i = numsize-ARGHEAD; NCOPY(t2,t1,i) }
2361 t1 = den1;
2362 if ( *den1 < 0 ) { *t2++ = den1[0]; *t2++ = den1[1]; }
2363 else { *t2++ = densize; *t2++ = 0; FILLARG(t2);
2364 i = densize-ARGHEAD; NCOPY(t2,t1,i) }
2365
2366 TermFree(num1,"MULfunc");
2367 TermFree(den1,"MULfunc");
2368 t = nextt;
2369 }
2370
2371 if ( flag ) AT.WorkPointer = term +*term;
2372 else AT.WorkPointer = oldworkpointer;
2373 AR.CompareRoutine = oldcompareroutine;
2374 AR.SortType = oldsorttype;
2375 return(0);
2376}
2377
2378/*
2379 #] ReadPolyRatFun :
2380 #[ FromPolyRatFun :
2381*/
2382
2383int 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);
2409Improper:
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*/
2433
2434WORD *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 }
2492didwork:;
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);
2623OverWork:
2624 MLOCK(ErrorMessageLock);
2625 MesPrint("Term too complex. Maybe increasing MaxTermSize can help");
2626 MUNLOCK(ErrorMessageLock);
2627CalledFrom:
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
2644void 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
2741WORD *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
2786WORD divrem[4] = { DIVFUNCTION, REMFUNCTION, INVERSEFUNCTION, MULFUNCTION };
2787
2788int DIVfunction(PHEAD WORD *term,WORD level,int par)
2789{
2790 GETBIDENTITY
2791 WORD *t, *tt, *r, *arg1 = 0, *arg2 = 0, *arg3 = 0, *termout;
2792 WORD *tstop, *tend, *r3, *rr, *rstop, tlength, rlength, newlength;
2793 WORD *proper1, *proper2, *proper3 = 0, numdol = -1;
2794 int numargs = 0, type1, type2, actionflag1, actionflag2;
2795 WORD startebuf = cbuf[AT.ebufnum].numrhs;
2796 int division = ( par <= 2 ); /* false for mul_ */
2797 if ( par < 0 || par > 3 ) {
2798 MLOCK(ErrorMessageLock);
2799 MesPrint("Internal error. Illegal parameter %d in DIVfunction.",par);
2800 MUNLOCK(ErrorMessageLock);
2801 Terminate(-1);
2802 }
2803/*
2804 Find the function
2805*/
2806 tend = term + *term; tstop = tend - ABS(tend[-1]);
2807 t = term+1;
2808 while ( t < tstop ) {
2809 if ( *t != divrem[par] ) { t += t[1]; continue; }
2810 r = t + FUNHEAD;
2811 tt = t + t[1]; numargs = 0;
2812 while ( r < tt ) {
2813 if ( numargs == 0 ) { arg1 = r; }
2814 if ( numargs == 1 ) { arg2 = r; }
2815 if ( numargs == 2 && *r == -DOLLAREXPRESSION ) { numdol = r[1]; }
2816 numargs++;
2817 NEXTARG(r);
2818 }
2819 if ( numargs == 2 ) break;
2820 if ( division && numargs == 3 ) break;
2821 t = tt;
2822 }
2823 if ( t >= tstop ) {
2824 MLOCK(ErrorMessageLock);
2825 MesPrint("Internal error. Indicated div_ or rem_ function not encountered.");
2826 MUNLOCK(ErrorMessageLock);
2827 Terminate(-1);
2828 }
2829/*
2830 We have two arguments in arg1 and arg2.
2831*/
2832 if ( division && *arg1 == -SNUMBER && arg1[1] == 0 ) {
2833 if ( *arg2 == -SNUMBER && arg2[1] == 0 ) {
2834zerozero:;
2835 MLOCK(ErrorMessageLock);
2836 MesPrint("0/0 in either div_ or rem_ function.");
2837 MUNLOCK(ErrorMessageLock);
2838 Terminate(-1);
2839 }
2840 if ( numdol >= 0 ) PutTermInDollar(0,numdol);
2841 return(0);
2842 }
2843 if ( division && *arg2 == -SNUMBER && arg2[1] == 0 ) {
2844divzero:;
2845 MLOCK(ErrorMessageLock);
2846 MesPrint("Division by zero in either div_ or rem_ function.");
2847 MUNLOCK(ErrorMessageLock);
2848 Terminate(-1);
2849 }
2850 if ( !division ) {
2851 if ( (*arg1 == -SNUMBER && arg1[1] == 0) ||
2852 (*arg2 == -SNUMBER && arg2[1] == 0) ) {
2853 return(0);
2854 }
2855 }
2856 if ( ( arg1 = ConvertArgument(BHEAD arg1, &type1) ) == 0 ) goto CalledFrom;
2857 if ( ( arg2 = ConvertArgument(BHEAD arg2, &type2) ) == 0 ) goto CalledFrom;
2858 if ( division && *arg1 == 0 ) {
2859 if ( *arg2 == 0 ) {
2860 M_free(arg2,"DIVfunction");
2861 M_free(arg1,"DIVfunction");
2862 goto zerozero;
2863 }
2864 M_free(arg2,"DIVfunction");
2865 M_free(arg1,"DIVfunction");
2866 if ( numdol >= 0 ) PutTermInDollar(0,numdol);
2867 return(0);
2868 }
2869 if ( division && *arg2 == 0 ) {
2870 M_free(arg2,"DIVfunction");
2871 M_free(arg1,"DIVfunction");
2872 goto divzero;
2873 }
2874 if ( !division && (*arg1 == 0 || *arg2 == 0) ) {
2875 M_free(arg2,"DIVfunction");
2876 M_free(arg1,"DIVfunction");
2877 return(0);
2878 }
2879 if ( ( proper1 = PutExtraSymbols(BHEAD arg1,startebuf,&actionflag1) ) == 0 ) goto CalledFrom;
2880 if ( ( proper2 = PutExtraSymbols(BHEAD arg2,startebuf,&actionflag2) ) == 0 ) goto CalledFrom;
2881/*
2882 if ( type2 == 0 ) M_free(arg2,"DIVfunction");
2883 else {
2884 DOLLARS d = ((DOLLARS)arg2)-1;
2885 if ( d->factors ) M_free(d->factors,"Dollar factors");
2886 M_free(d,"Copy of dollar variable");
2887 }
2888*/
2889 M_free(arg2,"DIVfunction");
2890/*
2891 if ( type1 == 0 ) M_free(arg1,"DIVfunction");
2892 else {
2893 DOLLARS d = ((DOLLARS)arg1)-1;
2894 if ( d->factors ) M_free(d->factors,"Dollar factors");
2895 M_free(d,"Copy of dollar variable");
2896 }
2897*/
2898 M_free(arg1,"DIVfunction");
2899 if ( par == 0 ) proper3 = poly_div(BHEAD proper1, proper2,0);
2900 else if ( par == 1 ) proper3 = poly_rem(BHEAD proper1, proper2,0);
2901 else if ( par == 2 ) proper3 = poly_inverse(BHEAD proper1, proper2);
2902 else if ( par == 3 ) proper3 = poly_mul(BHEAD proper1, proper2);
2903 if ( proper3 == 0 ) goto CalledFrom;
2904 if ( actionflag1 || actionflag2 ) {
2905 if ( ( arg3 = TakeExtraSymbols(BHEAD proper3,startebuf) ) == 0 ) goto CalledFrom;
2906 M_free(proper3,"DIVfunction");
2907 }
2908 else {
2909 arg3 = proper3;
2910 }
2911 M_free(proper2,"DIVfunction");
2912 M_free(proper1,"DIVfunction");
2913 cbuf[AT.ebufnum].numrhs = startebuf;
2914 if ( *arg3 ) {
2915 termout = AT.WorkPointer;
2916 tlength = tend[-1];
2917 tlength = REDLENG(tlength);
2918 r3 = arg3;
2919 while ( *r3 ) {
2920 tt = term + 1; rr = termout + 1;
2921 while ( tt < t ) *rr++ = *tt++;
2922 r = r3 + 1;
2923 r3 = r3 + *r3;
2924 rstop = r3 - ABS(r3[-1]);
2925 while ( r < rstop ) *rr++ = *r++;
2926 tt += t[1];
2927 while ( tt < tstop ) *rr++ = *tt++;
2928 rlength = r3[-1];
2929 rlength = REDLENG(rlength);
2930 if ( MulRat(BHEAD (UWORD *)tstop,tlength,(UWORD *)rstop,rlength,
2931 (UWORD *)rr,&newlength) < 0 ) goto CalledFrom;
2932 rlength = INCLENG(newlength);
2933 rr += ABS(rlength);
2934 rr[-1] = rlength;
2935 *termout = rr - termout;
2936 AT.WorkPointer = rr;
2937 if ( Generator(BHEAD termout,level) ) goto CalledFrom;
2938 }
2939 AT.WorkPointer = termout;
2940 }
2941 M_free(arg3,"DIVfunction");
2942 return(0);
2943CalledFrom:
2944 MLOCK(ErrorMessageLock);
2945 MesCall("DIVfunction");
2946 MUNLOCK(ErrorMessageLock);
2947 SETERROR(-1)
2948}
2949
2950/*
2951 #] DIVfunction :
2952 #[ MULfunc :
2953
2954 Multiplies two polynomials and puts the results in TermMalloc space.
2955*/
2956
2957WORD *MULfunc(PHEAD WORD *p1, WORD *p2)
2958{
2959 WORD *prod,size1,size2,size3,*t,*tfill,*ps1,*ps2,sign1,sign2, error, *p3;
2960 UWORD *num1, *num2, *num3;
2961 int i;
2962 WORD oldsorttype = AR.SortType;
2963 COMPARE oldcompareroutine = AR.CompareRoutine;
2964 AR.SortType = SORTHIGHFIRST;
2965 AR.CompareRoutine = &CompareSymbols;
2966 num3 = NumberMalloc("MULfunc");
2967 prod = TermMalloc("MULfunc");
2968 NewSort(BHEAD0);
2969 while ( *p1 ) {
2970 ps1 = p1+*p1; num1 = (UWORD *)(ps1 - ABS(ps1[-1])); size1 = ps1[-1];
2971 if ( size1 < 0 ) { sign1 = -1; size1 = -size1; }
2972 else sign1 = 1;
2973 size1 = (size1-1)/2;
2974 p3 = p2;
2975 while ( *p3 ) {
2976 ps2 = p3+*p3; num2 = (UWORD *)(ps2 - ABS(ps2[-1])); size2 = ps2[-1];
2977 if ( size2 < 0 ) { sign2 = -1; size2 = -size2; }
2978 else sign2 = 1;
2979 size2 = (size2-1)/2;
2980 if ( MulLong(num1,size1,num2,size2,num3,&size3) ) {
2981 error = 1;
2982CalledFrom:
2983 MLOCK(ErrorMessageLock);
2984 MesPrint(" Error %d",error);
2985 MesCall("MulFunc");
2986 MUNLOCK(ErrorMessageLock);
2987 Terminate(-1);
2988 }
2989 tfill = prod+1;
2990 t = p1+1; while ( t < (WORD *)num1 ) *tfill++ = *t++;
2991 t = p3+1; while ( t < (WORD *)num2 ) *tfill++ = *t++;
2992 t = (WORD *)num3;
2993 for ( i = 0; i < size3; i++ ) *tfill++ = *t++;
2994 *tfill++ = 1;
2995 for ( i = 1; i < size3; i++ ) *tfill++ = 0;
2996 *tfill++ = (2*size3+1)*sign1*sign2;
2997 prod[0] = tfill - prod;
2998 if ( SymbolNormalize(prod) ) { error = 2; goto CalledFrom; }
2999 if ( StoreTerm(BHEAD prod) ) { error = 3; goto CalledFrom; }
3000 p3 += *p3;
3001 }
3002 p1 += *p1;
3003 }
3004 NumberFree(num3,"MULfunc");
3005 EndSort(BHEAD prod,1);
3006 AR.CompareRoutine = oldcompareroutine;
3007 AR.SortType = oldsorttype;
3008 return(prod);
3009}
3010
3011/*
3012 #] MULfunc :
3013 #[ ConvertArgument :
3014
3015 Converts an argument to a general notation in allocated space.
3016*/
3017
3018WORD *ConvertArgument(PHEAD WORD *arg, int *type)
3019{
3020 WORD *output, *t, *r;
3021 int i, size;
3022 if ( *arg > 0 ) {
3023 output = (WORD *)Malloc1((*arg)*sizeof(WORD),"ConvertArgument");
3024 i = *arg - ARGHEAD; t = arg + ARGHEAD; r = output;
3025 NCOPY(r,t,i);
3026 *r = 0; *type = 0;
3027 return(output);
3028 }
3029 if ( *arg == -EXPRESSION ) {
3030 *type = 0;
3031 return(CreateExpression(BHEAD arg[1]));
3032 }
3033 if ( *arg == -DOLLAREXPRESSION ) {
3034 DOLLARS d;
3035 *type = 1;
3036 d = DolToTerms(BHEAD arg[1]);
3037/*
3038 The problem is that DolToTerms creates a copy of the dollar variable.
3039 If we just return d->where we create a memory leak. Hence we have to
3040 copy the contents of d->where to a new buffer
3041*/
3042 output = (WORD *)Malloc1((d->size+1)*sizeof(WORD),"Copy of dollar content");
3043 WCOPY(output,d->where,d->size+1);
3044 if ( d->factors ) { M_free(d->factors,"Dollar factors"); d->factors = 0; }
3045 M_free(d,"Copy of dollar variable");
3046 return(output);
3047 }
3048#if ( FUNHEAD > 4 )
3049 size = FUNHEAD+5;
3050#else
3051 size = 9;
3052#endif
3053 output = (WORD *)Malloc1(size*sizeof(WORD),"ConvertArgument");
3054 switch(*arg) {
3055 case -SYMBOL:
3056 output[0] = 8; output[1] = SYMBOL; output[2] = 4; output[3] = arg[1];
3057 output[4] = 1; output[5] = 1; output[6] = 1; output[7] = 3;
3058 output[8] = 0;
3059 break;
3060 case -INDEX:
3061 case -VECTOR:
3062 case -MINVECTOR:
3063 output[0] = 7; output[1] = INDEX; output[2] = 3; output[3] = arg[1];
3064 output[4] = 1; output[5] = 1;
3065 if ( *arg == -MINVECTOR ) output[6] = -3;
3066 else output[6] = 3;
3067 output[7] = 0;
3068 break;
3069 case -SNUMBER:
3070 output[0] = 4;
3071 if ( arg[1] < 0 ) {
3072 output[1] = -arg[1]; output[2] = 1; output[3] = -3;
3073 }
3074 else {
3075 output[1] = arg[1]; output[2] = 1; output[3] = 3;
3076 }
3077 output[4] = 0;
3078 break;
3079 default:
3080 output[0] = FUNHEAD+4;
3081 output[1] = -*arg;
3082 output[2] = FUNHEAD;
3083 for ( i = 3; i <= FUNHEAD; i++ ) output[i] = 0;
3084 output[FUNHEAD+1] = 1;
3085 output[FUNHEAD+2] = 1;
3086 output[FUNHEAD+3] = 3;
3087 output[FUNHEAD+4] = 0;
3088 break;
3089 }
3090 *type = 0;
3091 return(output);
3092}
3093
3094/*
3095 #] ConvertArgument :
3096 #[ ExpandRat :
3097
3098 Expands the denominator of a PolyRatFun in the variable PolyFunVar.
3099 The output is a polyratfun with a single argument.
3100 In the case that there is a polyratfun with more than one argument
3101 or the dirtyflag is on, the argument(s) is/are normalized.
3102 The output overwrites the input.
3103*/
3104
3105char *TheErrorMessage[] = {
3106 "PolyRatFun not of a type that FORM will expand: incorrect variable inside."
3107 ,"Division by zero in PolyRatFun encountered in ExpandRat."
3108 ,"Irregular code in PolyRatFun encountered in ExpandRat."
3109 ,"Called from ExpandRat."
3110 ,"WorkSpace overflow. Change parameter WorkSpace in setup file?"
3111 };
3112
3113int ExpandRat(PHEAD WORD *fun)
3114{
3115 WORD *r, *rr, *rrr, *tt, *tnext, *arg1, *arg2, *rmin = 0, *rmininv;
3116 WORD *rcoef, rsize, rcopy, *ow = AT.WorkPointer;
3117 WORD *numerator, *denominator, *rnext;
3118 WORD *thecopy, *rc, ncoef, newcoef, *m, *mm, nco, *outarg = 0;
3119 UWORD co[2], co1[2], co2[2];
3120 WORD OldPolyFunPow = AR.PolyFunPow;
3121 int i, j, minpow = 0, eppow, first, error = 0, ipoly;
3122 if ( fun[1] == FUNHEAD ) { return(0); }
3123 tnext = fun + fun[1];
3124 if ( fun[1] == fun[FUNHEAD]+FUNHEAD ) { /* Single argument */
3125 if ( fun[2] == 0 ) { goto done; }
3126/*
3127 We have to normalize the argument. This could make it shorter.
3128*/
3129NormArg:;
3130 if ( outarg == 0 ) outarg = TermMalloc("ExpandRat")+ARGHEAD;
3131 AT.TrimPower = 1;
3132 NewSort(BHEAD0);
3133 r = fun+FUNHEAD+ARGHEAD;
3134 if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3135 WORD minpow2 = MAXPOWER, *rrm;
3136 rrm = r;
3137 while ( rrm < tnext ) {
3138 if ( *rrm == 4 ) {
3139 if ( minpow2 > 0 ) minpow2 = 0;
3140 }
3141 else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3142 if ( minpow2 > 0 ) minpow2 = 0;
3143 }
3144 else {
3145 if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3146 if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3147 }
3148 else {
3149 MesPrint("Illegal term in expanded polyratfun.");
3150 goto onerror;
3151 }
3152 }
3153 rrm += *rrm;
3154 }
3155 AR.PolyFunPow += minpow2;
3156 }
3157 while ( r < tnext ) {
3158 rr = r + *r;
3159 i = *r; rrr = outarg; NCOPY(rrr,r,i);
3160 Normalize(BHEAD outarg);
3161 if ( *outarg > 0 ) StoreTerm(BHEAD outarg);
3162 }
3163 r = fun+FUNHEAD+ARGHEAD;
3164 EndSort(BHEAD r,1);
3165 AT.TrimPower = 0;
3166 if ( *r == 0 ) {
3167 fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0;
3168 fun[1] = FUNHEAD+2;
3169 }
3170 else {
3171 rr = fun+FUNHEAD;
3172 if ( ToFast(rr,rr) ) {
3173 NEXTARG(rr); fun[1] = rr - fun;
3174 }
3175 else {
3176 while ( *r ) r += *r;
3177 *rr = r-rr; rr[1] = CLEANFLAG;
3178 fun[1] = r - fun;
3179 }
3180 }
3181 fun[2] = CLEANFLAG;
3182 goto done;
3183 }
3184/*
3185 First test whether we have only AR.PolyFunVar in the denominator
3186*/
3187 tt = fun + FUNHEAD;
3188 arg1 = arg2 = 0;
3189 if ( tt < tnext ) {
3190 arg1 = tt; NEXTARG(tt);
3191 if ( tt < tnext ) {
3192 arg2 = tt; NEXTARG(tt);
3193 if ( tt != tnext ) { arg1 = arg2 = 0; } /* more than two arguments */
3194 }
3195 }
3196 if ( arg2 == 0 ) {
3197 if ( *arg1 < 0 ) { fun[2] = CLEANFLAG; goto done; }
3198 if ( fun[2] == CLEANFLAG ) goto done;
3199 goto NormArg; /* Note: should not come here */
3200 }
3201/*
3202 Produce the output argument in outarg
3203*/
3204 if ( outarg == 0 ) outarg = TermMalloc("ExpandRat")+ARGHEAD;
3205
3206 if ( *arg2 <= 0 ) {
3207/*
3208 These cases are trivial.
3209 We try as much as possible to write the output directly into the
3210 function. We just have to be extremely careful not to overwrite
3211 relevant information before we are finished with it.
3212*/
3213 if ( *arg2 == -SYMBOL && arg2[1] == AR.PolyFunVar ) {
3214 rr = r = fun+FUNHEAD+ARGHEAD;
3215 if ( *arg1 < 0 ) {
3216 if ( *arg1 == -SYMBOL ) {
3217 if ( arg1[1] == AR.PolyFunVar ) {
3218 *r++ = 4; *r++ = 1; *r++ = 1; *r++ = 3; *r = 0;
3219 }
3220 else {
3221 *r++ = 10; *r++ = SYMBOL; *r++ = 6;
3222 *r++ = arg1[1]; *r++ = 1;
3223 *r++ = AR.PolyFunVar; *r++ = -1;
3224 *r++ = 1; *r++ = 1; *r++ = 3; *r = 0;
3225 Normalize(BHEAD rr);
3226 }
3227 }
3228 else if ( *arg1 == -SNUMBER ) {
3229 nco = arg1[1];
3230 if ( nco == 0 ) { *r++ = 0; }
3231 else {
3232 *r++ = 8; *r++ = SYMBOL; *r++ = 4;
3233 *r++ = AR.PolyFunVar; *r++ = -1;
3234 *r++ = ABS(nco); *r++ = 1;
3235 if ( nco < 0 ) *r++ = -3;
3236 else *r++ = 3;
3237 *r = 0;
3238 }
3239 }
3240 else { error = 2; goto onerror; } /* should not happen! */
3241 }
3242 else { /* Multi-term numerator. */
3243 m = arg1+ARGHEAD;
3244 NewSort(BHEAD0); /* Technically maybe not needed */
3245 if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3246 WORD minpow2 = MAXPOWER, *rrm;
3247 rrm = m;
3248 while ( rrm < arg2 ) {
3249 if ( *rrm == 4 ) {
3250 if ( minpow2 > 0 ) minpow2 = 0;
3251 }
3252 else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3253 if ( minpow2 > 0 ) minpow2 = 0;
3254 }
3255 else {
3256 if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3257 if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3258 }
3259 else {
3260 MesPrint("Illegal term in expanded polyratfun.");
3261 goto onerror;
3262 }
3263 }
3264 rrm += *rrm;
3265 }
3266 AR.PolyFunPow += minpow2-1;
3267 }
3268 while ( m < arg2 ) {
3269 r = outarg;
3270 rrr = r++; mm = m + *m;
3271 *r++ = SYMBOL; *r++ = 4; *r++ = AR.PolyFunVar; *r++ = -1;
3272 m++; while ( m < mm ) *r++ = *m++;
3273 *rrr = r-rrr;
3274 Normalize(BHEAD rrr);
3275 StoreTerm(BHEAD rrr);
3276 }
3277 EndSort(BHEAD rr,1);
3278 r = rr; while ( *r ) r += *r;
3279 }
3280 if ( *rr == 0 ) {
3281 fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = CLEANFLAG;
3282 fun[1] = FUNHEAD+2;
3283 }
3284 else {
3285 rr = fun+FUNHEAD;
3286 *rr = r-rr;
3287 rr[1] = CLEANFLAG;
3288 if ( ToFast(rr,rr) ) {
3289 NEXTARG(rr);
3290 fun[1] = rr - fun;
3291 }
3292 else { fun[1] = r - fun; }
3293 }
3294 fun[2] = CLEANFLAG;
3295 goto done;
3296 }
3297 else if ( *arg2 == -SNUMBER ) {
3298 rr = r = outarg;
3299 if ( arg2[1] == 0 ) { error = 1; goto onerror; }
3300 if ( *arg1 == -SNUMBER ) { /* Things may not be normalized */
3301 if ( arg1[1] == 0 ) { *r++ = 0; }
3302 else {
3303 co1[0] = ABS(arg1[1]); co1[1] = 1;
3304 co2[0] = 1; co2[1] = ABS(arg2[1]);
3305 MulRat(BHEAD co1,1,co2,1,co,&nco);
3306 *r++ = 4; *r++ = (WORD)(co[0]); *r++ = (WORD)(co[1]);
3307 if ( ( arg1[1] < 0 && arg2[1] > 0 ) ||
3308 ( arg1[1] > 0 && arg2[1] < 0 ) ) *r++ = -3;
3309 else *r++ = 3;
3310 *r = 0;
3311 }
3312 }
3313 else if ( *arg1 == -SYMBOL ) {
3314 *r++ = 8; *r++ = SYMBOL; *r++ = 4;
3315 *r++ = arg1[1]; *r++ = 1;
3316 *r++ = 1; *r++ = ABS(arg2[1]);
3317 if ( arg2[1] < 0 ) *r++ = -3;
3318 else *r++ = 3;
3319 *r = 0;
3320 }
3321 else if ( *arg1 < 0 ) { error = 2; goto onerror; }
3322 else { /* Multi-term numerator. */
3323 m = arg1+ARGHEAD;
3324 NewSort(BHEAD0); /* Technically maybe not needed */
3325 if ( AR.PolyFunExp == 2 ) { /* Find minimum power */
3326 WORD minpow2 = MAXPOWER, *rrm;
3327 rrm = m;
3328 while ( rrm < arg2 ) {
3329 if ( *rrm == 4 ) {
3330 if ( minpow2 > 0 ) minpow2 = 0;
3331 }
3332 else if ( ABS(rrm[*rrm-1]) == (*rrm-1) ) {
3333 if ( minpow2 > 0 ) minpow2 = 0;
3334 }
3335 else {
3336 if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) {
3337 if ( rrm[4] < minpow2 ) minpow2 = rrm[4];
3338 }
3339 else {
3340 MesPrint("Illegal term in expanded polyratfun.");
3341 goto onerror;
3342 }
3343 }
3344 rrm += *rrm;
3345 }
3346 AR.PolyFunPow += minpow2;
3347 }
3348 while ( m < arg2 ) {
3349 r = rr;
3350 rrr = r++; mm = m + *m;
3351 *r++ = DENOMINATOR; *r++ = FUNHEAD + 2; *r++ = DIRTYFLAG;
3352 FILLFUN3(r);
3353 *r++ = arg2[0]; *r++ = arg2[1];
3354 m++; while ( m < mm ) *r++ = *m++;
3355 *rrr = r-rrr;
3356 if ( r < AT.WorkTop && r >= AT.WorkSpace )
3357 AT.WorkPointer = r;
3358 Normalize(BHEAD rrr);
3359 if ( ABS(rrr[*rrr-1]) == *rrr-1 ) {
3360 if ( AR.PolyFunPow >= 0 ) {
3361 StoreTerm(BHEAD rrr);
3362 }
3363 }
3364 else if ( rrr[1] == SYMBOL && rrr[2] == 4 &&
3365 rrr[3] == AR.PolyFunVar && rrr[4] <= AR.PolyFunPow ) {
3366 StoreTerm(BHEAD rrr);
3367 }
3368 }
3369 EndSort(BHEAD rr,1);
3370 }
3371 r = rr; while ( *r ) r += *r;
3372 i = r-rr;
3373 r = fun + FUNHEAD + ARGHEAD;
3374 NCOPY(r,rr,i);
3375 rr = fun + FUNHEAD;
3376 *rr = r - rr; rr[1] = CLEANFLAG;
3377 if ( ToFast(rr,rr) ) {
3378 NEXTARG(rr);
3379 fun[1] = rr - fun;
3380 }
3381 else { fun[1] = r - fun; }
3382 fun[2] = CLEANFLAG;
3383 goto done;
3384 }
3385 else { error = 0; goto onerror; }
3386 }
3387 else {
3388 r = arg2+ARGHEAD; /* The argument ends at tnext */
3389 first = 1;
3390 while ( r < tnext ) {
3391 rr = r + *r; rr -= ABS(rr[-1]);
3392 if ( r+1 == rr ) {
3393 if ( first ) { minpow = 0; first = 0; rmin = r; }
3394 else if ( minpow > 0 ) { minpow = 0; rmin = r; }
3395 }
3396 else if ( r[1] != SYMBOL || r[2] != 4 || r[3] != AR.PolyFunVar
3397 || r[4] > MAXPOWER ) { error = 0; goto onerror; }
3398 else if ( first ) { minpow = r[4]; first = 0; rmin = r; }
3399 else if ( r[4] < minpow ) { minpow = r[4]; rmin = r; }
3400 r += *r;
3401 }
3402/*
3403 We have now:
3404 1: a numerator in arg1 which can contain several variables.
3405 2: a denominator in arg2 with at most only AR.PolyFunVar (ep).
3406 3: the minimum power in the denominator is minpow and the
3407 term with that minimum power is in rmin.
3408 Divide numerator and denominator by this minimum power.
3409 Determine the power range in the numerator.
3410 Call InvPoly.
3411 Multiply by the inverse in such a way that we never take more
3412 powers of ep than necessary.
3413*/
3414/*
3415 One: put 1/rmin in AT.WorkPointer -> rmininv
3416*/
3417 AT.WorkPointer += AM.MaxTer/sizeof(WORD);
3418 if ( AT.WorkPointer + (AM.MaxTer/sizeof(WORD)) >= AT.WorkTop ) {
3419 error = 4; goto onerror;
3420 }
3421 rmininv = r = AT.WorkPointer;
3422 rr = rmin; i = *rmin; NCOPY(r,rr,i)
3423 if ( minpow != 0 ) { rmininv[4] = -rmininv[4]; }
3424 rsize = ABS(r[-1]);
3425 rcoef = r - rsize;
3426 rsize = (rsize-1)/2; rr = rcoef + rsize;
3427 for ( i = 0; i < rsize; i++ ) {
3428 rcopy = rcoef[i]; rcoef[i] = rr[i]; rr[i] = rcopy;
3429 }
3430 AT.WorkPointer = r;
3431 if ( *arg1 < 0 ) {
3432 ToGeneral(arg1,r,0);
3433 arg1 = r; r += *r; *r++ = 0; rcopy = 0;
3434 AT.WorkPointer = r;
3435 }
3436 else {
3437 r = arg1 + *arg1;
3438 rcopy = *r; *r++ = 0;
3439 }
3440/*
3441 We can use MultiplyWithTerm.
3442*/
3443 AT.LeaveNegative = 1;
3444 numerator = MultiplyWithTerm(BHEAD arg1+ARGHEAD,rmininv,0);
3445 AT.LeaveNegative = 0;
3446 r[-1] = rcopy;
3447 r = numerator; while ( *r ) r += *r;
3448 AT.WorkPointer = r+1;
3449 rcopy = arg2[*arg2]; arg2[*arg2] = 0;
3450 denominator = MultiplyWithTerm(BHEAD arg2+ARGHEAD,rmininv,1);
3451 arg2[*arg2] = rcopy;
3452 r = denominator; while ( *r ) r += *r;
3453 AT.WorkPointer = r+1;
3454/*
3455 Now find the minimum power of ep in the numerator.
3456*/
3457 r = numerator;
3458 first = 1;
3459 while ( *r ) {
3460 rr = r + *r; rr -= ABS(rr[-1]);
3461 if ( r+1 == rr ) {
3462 if ( first ) { minpow = 0; first = 0; }
3463 else if ( minpow > 0 ) { minpow = 0; }
3464 }
3465 else if ( r[1] != SYMBOL ) { error = 0; goto onerror; }
3466 else {
3467 for ( i = 3; i < r[2]; i += 2 ) {
3468 if ( r[i] == AR.PolyFunVar ) {
3469 if ( first ) { minpow = r[i+1]; first = 0; }
3470 else if ( r[i+1] < minpow ) minpow = r[i+1];
3471 break;
3472 }
3473 }
3474 if ( i >= r[2] ) {
3475 if ( first ) { minpow = 0; first = 0; }
3476 else if ( minpow > 0 ) minpow = 0;
3477 }
3478 }
3479 r += *r;
3480 }
3481/*
3482 We can invert the denominator.
3483 Note that the return value is an offset in AT.pWorkSpace.
3484 Hence there is no need to free memory afterwards.
3485*/
3486 if ( AR.PolyFunExp == 3 ) {
3487 ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow-minpow,AR.PolyFunVar);
3488 }
3489 else {
3490 ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow,AR.PolyFunVar);
3491 }
3492/*
3493 Now we start the multiplying
3494*/
3495 NewSort(BHEAD0);
3496 r = numerator;
3497 while ( *r ) {
3498/*
3499 1: Find power of ep.
3500*/
3501 rnext = r + *r;
3502 rrr = rnext - ABS(rnext[-1]);
3503 rr = r+1;
3504 eppow = 0;
3505 if ( rr < rrr ) {
3506 j = rr[1] - 2; rr += 2;
3507 while ( j > 0 ) {
3508 if ( *rr == AR.PolyFunVar ) { eppow = rr[1]; break; }
3509 j -= 2; rr += 2;
3510 }
3511 }
3512/*
3513 2: Multiply by the proper terms in ipoly
3514*/
3515 for ( i = 0; i <= AR.PolyFunPow-eppow+minpow; i++ ) {
3516 if ( AT.pWorkSpace[ipoly+i] == 0 ) continue;
3517/*
3518 Copy the term, add i to the power of ep and multiply coef.
3519*/
3520 rc = r;
3521 rr = thecopy = AT.WorkPointer;
3522 while ( rc < rrr ) *rr++ = *rc++;
3523 if ( i != 0 ) {
3524 *rr++ = SYMBOL; *rr++ = 4; *rr++ = AR.PolyFunVar; *rr++ = i;
3525 }
3526 ncoef = REDLENG(rnext[-1]);
3527 MulRat(BHEAD (UWORD *)rrr,ncoef,
3528 (UWORD *)(AT.pWorkSpace[ipoly+i])+1,AT.pWorkSpace[ipoly+i][0]
3529 ,(UWORD *)rr,&newcoef);
3530 ncoef = ABS(newcoef); rr += 2*ncoef;
3531 newcoef = INCLENG(newcoef);
3532 *rr++ = newcoef;
3533 *thecopy = rr - thecopy;
3534 AT.WorkPointer = rr;
3535 Normalize(BHEAD thecopy);
3536 if ( *thecopy > 0 ) StoreTerm(BHEAD thecopy);
3537 AT.WorkPointer = thecopy;
3538 }
3539 r = rnext;
3540 }
3541/*
3542 Now we have all.
3543*/
3544 rr = fun + FUNHEAD; r = rr + ARGHEAD;
3545 EndSort(BHEAD r,1);
3546 if ( *r == 0 ) {
3547 fun[1] = FUNHEAD+2; fun[2] = CLEANFLAG;
3548 fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0;
3549 }
3550 else {
3551 while ( *r ) r += *r;
3552 rr[0] = r-rr; rr[1] = CLEANFLAG;
3553 if ( ToFast(rr,rr) ) { NEXTARG(rr); fun[1] = rr-fun; }
3554 else { fun[1] = r-fun; }
3555 fun[2] = CLEANFLAG;
3556 }
3557 }
3558done:
3559 if ( outarg ) TermFree(outarg-ARGHEAD,"ExpandRat");
3560 AR.PolyFunPow = OldPolyFunPow;
3561 AT.WorkPointer = ow;
3562 AN.PolyNormFlag = 1;
3563 return(0);
3564onerror:
3565 if ( outarg ) TermFree(outarg-ARGHEAD,"ExpandRat");
3566 AR.PolyFunPow = OldPolyFunPow;
3567 AT.WorkPointer = ow;
3568 MLOCK(ErrorMessageLock);
3569 MesPrint(TheErrorMessage[error]);
3570 MUNLOCK(ErrorMessageLock);
3571 Terminate(-1);
3572 return(-1);
3573}
3574
3575/*
3576 #] ExpandRat :
3577 #[ InvPoly :
3578
3579 The input polynomial is represented as a sequence of terms in ascending
3580 power. The first coefficient is 1. If we call this 1-a and
3581 a = sum_(j,1,n,x^j*a(j)), and b = 1/(1-a) we can find the coefficients
3582 of b with the recursion
3583 b(0) = 1, b(n) = sum_(j,1,n,a(j)*b(n-j))
3584 The variable is the symbol sym and we need maxpow powers in the answer.
3585 The answer is an array of pointers to the coefficients of the various
3586 powers as rational numbers in the notation signedsize,numerator,denominator
3587 We put these powers in the workspace and the answer is in AT.pWorkSpace.
3588 Hence the return value is an offset in the pWorkSpace.
3589 A zero pointer indicates that this coefficient is zero.
3590*/
3591
3592int InvPoly(PHEAD WORD *inpoly, WORD maxpow, WORD sym)
3593{
3594 int needed, inpointers, outpointers, maxinput = 0, i, j;
3595 WORD *t, *tt, *ttt, *w, *c, *cc, *ccc, lenc, lenc1, lenc2, rc, *c1, *c2;
3596/*
3597 Step 0: allocate the space
3598*/
3599 needed = (maxpow+1)*2;
3600 WantAddPointers(needed);
3601 inpointers = AT.pWorkPointer;
3602 outpointers = AT.pWorkPointer+maxpow+1;
3603 for ( i = 0; i < needed; i++ ) AT.pWorkSpace[inpointers+i] = 0;
3604/*
3605 Step 1: determine the coefficients in inpoly
3606 often there is a maximum power that is much smaller than maxpow.
3607 keeping track of this can speed up things.
3608*/
3609 t = inpoly;
3610 w = AT.WorkPointer;
3611 while ( *t ) {
3612 if ( *t == 4 ) {
3613 if ( t[1] != 1 || t[2] != 1 || t[3] != 3 ) goto onerror;
3614 AT.pWorkSpace[inpointers] = 0;
3615 }
3616 else if ( t[1] != SYMBOL || t[2] != 4 || t[3] != sym || t[4] < 0 ) goto onerror;
3617 else if ( t[4] > maxpow ) {} /* power outside useful range */
3618 else {
3619 if ( t[4] > maxinput ) maxinput = t[4];
3620 AT.pWorkSpace[inpointers+t[4]] = w;
3621 tt = t + *t; rc = -*--tt; /* we need - the coefficient! */
3622 rc = REDLENG(rc); *w++ = rc;
3623 ttt = t+5;
3624 while ( ttt < tt ) *w++ = *ttt++;
3625 }
3626 t += *t;
3627 }
3628/*
3629 Step 2: compute the output. b(0) = 1.
3630 then the recursion starts.
3631*/
3632 AT.pWorkSpace[outpointers] = w;
3633 *w++ = 1; *w++ = 1; *w++ = 1;
3634 c = TermMalloc("InvPoly");
3635 c1 = TermMalloc("InvPoly");
3636 c2 = TermMalloc("InvPoly");
3637 for ( j = 1; j <= maxpow; j++ ) {
3638/*
3639 Start at c = a(j)*b(0) = a(j)
3640*/
3641 if ( ( cc = AT.pWorkSpace[inpointers+j] ) != 0 ) {
3642 lenc = *cc++; /* reduced length */
3643 i = 2*ABS(lenc); ccc = c;
3644 NCOPY(ccc,cc,i);
3645 }
3646 else { lenc = 0; }
3647 for ( i = MiN(j-1,maxinput); i > 0; i-- ) {
3648/*
3649 c -> c + a(i)*b(j-i)
3650*/
3651 if ( AT.pWorkSpace[inpointers+i] == 0
3652 || AT.pWorkSpace[outpointers+j-i] == 0 ) {
3653 }
3654 else {
3655 if ( MulRat(BHEAD (UWORD *)(AT.pWorkSpace[inpointers+i]+1),AT.pWorkSpace[inpointers+i][0],
3656 (UWORD *)(AT.pWorkSpace[outpointers+j-i]+1),AT.pWorkSpace[outpointers+j-i][0],
3657 (UWORD *)c1,&lenc1) ) goto calcerror;
3658 if ( lenc == 0 ) {
3659 cc = c; c = c1; c1 = cc;
3660 lenc = lenc1;
3661 }
3662 else {
3663 if ( AddRat(BHEAD (UWORD *)c,lenc,(UWORD *)c1,lenc1,(UWORD *)c2,&lenc2) )
3664 goto calcerror;
3665 cc = c; c = c2; c2 = cc;
3666 lenc = lenc2;
3667 }
3668 }
3669 }
3670/*
3671 Copy c to the proper location
3672*/
3673 if ( lenc == 0 ) AT.pWorkSpace[outpointers+j] = 0;
3674 else {
3675 AT.pWorkSpace[outpointers+j] = w;
3676 *w++ = lenc;
3677 i = 2*ABS(lenc); ccc = c;
3678 NCOPY(w,ccc,i);
3679 }
3680 }
3681 AT.WorkPointer = w;
3682 TermFree(c2,"InvPoly");
3683 TermFree(c1,"InvPoly");
3684 TermFree(c ,"InvPoly");
3685
3686 return(outpointers);
3687onerror:
3688 MLOCK(ErrorMessageLock);
3689 MesPrint("Incorrect symbol field in InvPoly.");
3690 MUNLOCK(ErrorMessageLock);
3691 Terminate(-1);
3692 return(-1);
3693calcerror:
3694 MLOCK(ErrorMessageLock);
3695 MesPrint("Called from InvPoly.");
3696 MUNLOCK(ErrorMessageLock);
3697 Terminate(-1);
3698 return(-1);
3699}
3700
3701/*
3702 #] InvPoly :
3703*/
WORD CompareSymbols(WORD *, WORD *, WORD)
Definition sort.c:2976
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
Definition notation.c:510
WORD NewSort(PHEAD0)
Definition sort.c:592
LONG EndSort(PHEAD WORD *, int)
Definition sort.c:682
WORD InsertTerm(PHEAD WORD *, WORD, WORD, WORD *, WORD *, WORD)
Definition proces.c:2579
WORD Generator(PHEAD WORD *, WORD)
Definition proces.c:3101
WORD StoreTerm(PHEAD WORD *)
Definition sort.c:4333
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition reken.c:1466
int SymbolNormalize(WORD *)
Definition normal.c:5014
WORD * poly_gcd(PHEAD WORD *, WORD *, WORD)
Definition polywrap.cc:124
WORD * TakeContent(PHEAD WORD *in, WORD *term)
Definition ratio.c:1376
WORD * TakeSymbolContent(PHEAD WORD *in, WORD *term)
Definition ratio.c:2434
VOID LowerSortLevel()
Definition sort.c:4727
WORD ** rhs
Definition structs.h:943
WORD * Buffer
Definition structs.h:939
struct CbUf CBUF
struct FiLe FILEHANDLE