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