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