FORM v5.0.0-35-g6318119
argument.c
Go to the documentation of this file.
1
7/* #[ License : */
8/*
9 * Copyright (C) 1984-2026 J.A.M. Vermaseren
10 * When using this file you are requested to refer to the publication
11 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
12 * This is considered a matter of courtesy as the development was paid
13 * for by FOM the Dutch physics granting agency and we would like to
14 * be able to track its scientific use to convince FOM of its value
15 * for the community.
16 *
17 * This file is part of FORM.
18 *
19 * FORM is free software: you can redistribute it and/or modify it under the
20 * terms of the GNU General Public License as published by the Free Software
21 * Foundation, either version 3 of the License, or (at your option) any later
22 * version.
23 *
24 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
25 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
26 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
27 * details.
28 *
29 * You should have received a copy of the GNU General Public License along
30 * with FORM. If not, see <http://www.gnu.org/licenses/>.
31 */
32/* #] License : */
33
34/*
35 #[ include : argument.c
36*/
37
38#include "form3.h"
39
40/*
41 #] include :
42 #[ execarg :
43
44 Executes the subset of statements in an argument environment.
45 The calling routine should be of the type
46 if ( C->lhs[level][0] == TYPEARG ) {
47 if ( execarg(term,level) ) goto GenCall;
48 level = C->lhs[level][2];
49 goto SkipCount;
50 }
51 Note that there will be cases in which extra space is needed.
52 In addition the compare with C->numlhs isn't very fine, because we
53 need to insert a different value (C->lhs[level][2]).
54*/
55
56int execarg(PHEAD WORD *term, WORD level)
57{
58 GETBIDENTITY
59 WORD *t, *r, *m, *v;
60 WORD *start, *stop, *rstop, *r1, *r2 = 0, *r3 = 0, *r4, *r5, *r6, *r7, *r8, *r9;
61 WORD *mm, *mstop, *rnext, *rr, *factor, type, ngcd, nq;
62 CBUF *C = cbuf+AM.rbufnum, *CC = cbuf+AT.ebufnum;
63 WORD i, j, k, oldnumlhs = AR.Cnumlhs, count, olddefer = AR.DeferFlag;
64 int action = 0;
65 WORD oldnumrhs = CC->numrhs, size, pow, jj;
66 LONG oldcpointer = CC->Pointer - CC->Buffer, oldppointer = AT.pWorkPointer, lp;
67 WORD *oldwork = AT.WorkPointer, *oldwork2, scale, renorm;
68 WORD kLCM = 0, kGCD = 0, kGCD2, kkLCM = 0, jLCM = 0, jGCD, sign = 1;
69 int ii, didpolyratfun;
70 UWORD *EAscrat, *GCDbuffer = 0, *GCDbuffer2 = 0, *LCMbuffer = 0, *LCMb = 0, *LCMc = 0;
71 AT.WorkPointer += *term;
72 start = C->lhs[level];
73 AR.Cnumlhs = start[2];
74 stop = start + start[1];
75 type = *start;
76 scale = start[4];
77 renorm = start[5];
78 start += TYPEARGHEADSIZE;
79/*
80 #[ Dollars :
81*/
82 if ( renorm && start[1] != 0 ) {/* We have to evaluate $ symbols inside () */
83 t = start+1; factor = oldwork2 = v = AT.WorkPointer;
84 i = *t; t++;
85 *v++ = i+3; i--; NCOPY(v,t,i);
86 *v++ = 1; *v++ = 1; *v++ = 3;
87 AT.WorkPointer = v;
88 start = t; AR.Eside = LHSIDEX;
89 NewSort(BHEAD0);
90 if ( Generator(BHEAD factor,AR.Cnumlhs) ) {
91 LowerSortLevel();
92 AT.WorkPointer = oldwork;
93 return(-1);
94 }
95 AT.WorkPointer = v;
96 if ( EndSort(BHEAD factor,0) < 0 ) {}
97 if ( *factor && *(factor+*factor) != 0 ) {
98 MLOCK(ErrorMessageLock);
99 MesPrint("&$ in () does not evaluate into a single term");
100 MUNLOCK(ErrorMessageLock);
101 return(-1);
102 }
103 AR.Eside = RHSIDE;
104 if ( *factor > 0 ) {
105 v = factor+*factor;
106 v -= ABS(v[-1]);
107 *factor = v-factor;
108 }
109 AT.WorkPointer = v;
110 }
111 else {
112 if ( *start < 0 ) {
113 factor = start + 1;
114 start += -*start;
115 }
116 else factor = 0;
117 }
118/*
119 #] Dollars :
120*/
121 t = term;
122 r = t + *t;
123 rstop = r - ABS(r[-1]);
124 t++;
125/*
126 #[ Argument detection : + argument statement
127*/
128/*
129 Allocate Numbers for MakeInteger here, for re-use in case multiple
130 functions are treated in the same term.
131*/
132 if ( type == TYPENORM4 ) {
133 GCDbuffer = NumberMalloc("execarg");
134 GCDbuffer2 = NumberMalloc("execarg");
135 LCMbuffer = NumberMalloc("execarg");
136 LCMb = NumberMalloc("execarg"); LCMc = NumberMalloc("execarg");
137 }
138 didpolyratfun = 0;
139 while ( t < rstop ) {
140 if ( *t >= FUNCTION && functions[*t-FUNCTION].spec <= 0 ) {
141/*
142 We have a function. First count the number of arguments.
143 Tensors are excluded.
144*/
145 count = 0;
146 v = t;
147 m = t + FUNHEAD;
148 r = t + t[1];
149 while ( m < r ) {
150 count++;
151 NEXTARG(m)
152 }
153 if ( count <= 0 ) { t += t[1]; continue; }
154/*
155 Now we take the arguments one by one and test for a match
156*/
157 for ( i = 1; i <= count; i++ ) {
158 m = start;
159 while ( m < stop ) {
160 r = m + m[1];
161 j = *r++;
162 if ( j > 1 ) {
163 while ( --j > 0 ) {
164 if ( *r == i ) goto RightNum;
165 r++;
166 }
167 m = r;
168 continue;
169 }
170RightNum:
171 if ( m[1] == 2 ) {
172#ifdef WITHFLOAT
173 if ( *t != FLOATFUN || TestFloat(t) == 0 )
174#endif
175 {
176 m += 2;
177 m += *m;
178 goto HaveTodo;
179 }
180#ifdef WITHFLOAT
181 else {
182 m += 2;
183 }
184#endif
185 }
186 else {
187 r = m + m[1];
188 m += 2;
189 while ( m < r ) {
190 if ( *m == CSET ) {
191 r1 = SetElements + Sets[m[1]].first;
192 r2 = SetElements + Sets[m[1]].last;
193 while ( r1 < r2 ) {
194 if ( *r1++ == *t ) goto HaveTodo;
195 }
196 }
197 else if ( m[1] == *t ) goto HaveTodo;
198 m += 2;
199 }
200 }
201 m += *m;
202 }
203 continue;
204HaveTodo:
205/*
206 If we come here we have to do the argument i (first is 1).
207*/
208 sign = 1;
209 action = 1;
210 if ( *t == AR.PolyFun ) didpolyratfun = 1;
211 v[2] |= DIRTYFLAG;
212 r = t + FUNHEAD;
213 j = i;
214 while ( --j > 0 ) { NEXTARG(r) }
215 if ( ( type == TYPESPLITARG ) || ( type == TYPESPLITFIRSTARG )
216 || ( type == TYPESPLITLASTARG ) ) {
217 if ( *t > FUNCTION && *r > 0 ) {
218 WantAddPointers(2);
219 AT.pWorkSpace[AT.pWorkPointer++] = t;
220 AT.pWorkSpace[AT.pWorkPointer++] = r;
221 }
222 continue;
223 }
224 else if ( type == TYPESPLITARG2 ) {
225 if ( *t > FUNCTION && *r > 0 ) {
226 WantAddPointers(2);
227 AT.pWorkSpace[AT.pWorkPointer++] = t;
228 AT.pWorkSpace[AT.pWorkPointer++] = r;
229 }
230 continue;
231 }
232 else if ( type == TYPEFACTARG || type == TYPEFACTARG2 ) {
233 if ( *t > FUNCTION || *t == DENOMINATOR ) {
234 if ( *r > 0 ) {
235 mm = r + ARGHEAD; mstop = r + *r;
236 if ( mm + *mm < mstop ) {
237 WantAddPointers(2);
238 AT.pWorkSpace[AT.pWorkPointer++] = t;
239 AT.pWorkSpace[AT.pWorkPointer++] = r;
240 continue;
241 }
242 if ( *mm == 1+ABS(mstop[-1]) ) continue;
243 if ( mstop[-3] != 1 || mstop[-2] != 1
244 || mstop[-1] != 3 ) {
245 WantAddPointers(2);
246 AT.pWorkSpace[AT.pWorkPointer++] = t;
247 AT.pWorkSpace[AT.pWorkPointer++] = r;
248 continue;
249 }
250 GETSTOP(mm,mstop); mm++;
251 if ( mm + mm[1] < mstop ) {
252 WantAddPointers(2);
253 AT.pWorkSpace[AT.pWorkPointer++] = t;
254 AT.pWorkSpace[AT.pWorkPointer++] = r;
255 continue;
256 }
257 if ( *mm == SYMBOL && ( mm[1] > 4 ||
258 ( mm[3] != 1 && mm[3] != -1 ) ) ) {
259 WantAddPointers(2);
260 AT.pWorkSpace[AT.pWorkPointer++] = t;
261 AT.pWorkSpace[AT.pWorkPointer++] = r;
262 continue;
263 }
264 else if ( *mm == DOTPRODUCT && ( mm[1] > 5 ||
265 ( mm[4] != 1 && mm[4] != -1 ) ) ) {
266 WantAddPointers(2);
267 AT.pWorkSpace[AT.pWorkPointer++] = t;
268 AT.pWorkSpace[AT.pWorkPointer++] = r;
269 continue;
270 }
271 else if ( ( *mm == DELTA || *mm == VECTOR )
272 && mm[1] > 4 ) {
273 WantAddPointers(2);
274 AT.pWorkSpace[AT.pWorkPointer++] = t;
275 AT.pWorkSpace[AT.pWorkPointer++] = r;
276 continue;
277 }
278 }
279 else if ( factor && *factor == 4 && factor[2] == 1 ) {
280 WantAddPointers(2);
281 AT.pWorkSpace[AT.pWorkPointer++] = t;
282 AT.pWorkSpace[AT.pWorkPointer++] = r;
283 continue;
284 }
285 else if ( factor && *factor == 0
286 && ( *r == -SNUMBER && r[1] != 1 ) ) {
287 WantAddPointers(2);
288 AT.pWorkSpace[AT.pWorkPointer++] = t;
289 AT.pWorkSpace[AT.pWorkPointer++] = r;
290 continue;
291 }
292 else if ( *r == -MINVECTOR ) {
293 WantAddPointers(2);
294 AT.pWorkSpace[AT.pWorkPointer++] = t;
295 AT.pWorkSpace[AT.pWorkPointer++] = r;
296 continue;
297 }
298 }
299 continue;
300 }
301 else if ( type == TYPENORM || type == TYPENORM2 || type == TYPENORM3 || type == TYPENORM4 ) {
302 if ( *r < 0 ) {
303 WORD rone;
304 if ( *r == -MINVECTOR ) { rone = -1; *r = -INDEX; }
305 else if ( *r != -SNUMBER || r[1] == 1 || r[1] == 0 ) continue;
306 else { rone = r[1]; r[1] = 1; }
307/*
308 Now we must multiply the general coefficient by r[1]
309*/
310 if ( scale && ( factor == 0 || *factor ) ) {
311 action = 1;
312 v[2] |= DIRTYFLAG;
313 if ( rone < 0 ) {
314 if ( type == TYPENORM3 ) k = 1;
315 else k = -1;
316 rone = -rone;
317 }
318 else k = 1;
319 r1 = term + *term;
320 size = r1[-1];
321 size = REDLENG(size);
322 if ( scale > 0 ) {
323 for ( jj = 0; jj < scale; jj++ ) {
324 if ( Mully(BHEAD (UWORD *)rstop,&size,(UWORD *)(&rone),k) )
325 goto execargerr;
326 }
327 }
328 else {
329 for ( jj = 0; jj > scale; jj-- ) {
330 if ( Divvy(BHEAD (UWORD *)rstop,&size,(UWORD *)(&rone),k) )
331 goto execargerr;
332 }
333 }
334 size = INCLENG(size);
335 k = size < 0 ? -size: size;
336 rstop[k-1] = size;
337 *term = (WORD)(rstop - term) + k;
338 }
339 continue;
340 }
341/*
342 Now we have to find a reference term.
343 If factor is defined and *factor != 0 we have to
344 look for the first term that matches the pattern exactly
345 Otherwise the first term plays this role
346 If its coefficient is not one,
347 we must set up a division of the whole argument by
348 this coefficient, and a multiplication of the term
349 when the type is not equal to TYPENORM2.
350 We first multiply the coefficient of the term.
351 Then we set up the division.
352
353 First find the magic term
354*/
355 if ( type == TYPENORM4 ) {
356/*
357 For normalizing everything to integers we have to
358 determine for all elements of this argument the LCM of
359 the denominators and the GCD of the numerators.
360 The buffers have been allocated already.
361*/
362 r4 = r + *r;
363 r1 = r + ARGHEAD;
364/*
365 First take the first term to load up the LCM and the GCD
366*/
367 r2 = r1 + *r1;
368 j = r2[-1];
369 if ( j < 0 ) sign = -1;
370 r3 = r2 - ABS(j);
371 k = REDLENG(j);
372 if ( k < 0 ) k = -k;
373 while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
374 for ( kGCD = 0; kGCD < k; kGCD++ ) GCDbuffer[kGCD] = r3[kGCD];
375 k = REDLENG(j);
376 if ( k < 0 ) k = -k;
377 r3 += k;
378 while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
379 for ( kLCM = 0; kLCM < k; kLCM++ ) LCMbuffer[kLCM] = r3[kLCM];
380 r1 = r2;
381/*
382 Now go through the rest of the terms in this argument.
383*/
384 while ( r1 < r4 ) {
385 r2 = r1 + *r1;
386 j = r2[-1];
387 r3 = r2 - ABS(j);
388 k = REDLENG(j);
389 if ( k < 0 ) k = -k;
390 while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
391 if ( ( ( GCDbuffer[0] == 1 ) && ( kGCD == 1 ) ) ) {
392/*
393 GCD is already 1
394*/
395 }
396 else if ( ( ( k != 1 ) || ( r3[0] != 1 ) ) ) {
397 if ( GcdLong(BHEAD GCDbuffer,kGCD,(UWORD *)r3,k,GCDbuffer2,&kGCD2) ) {
398 NumberFree(GCDbuffer,"execarg");
399 NumberFree(GCDbuffer2,"execarg");
400 NumberFree(LCMbuffer,"execarg");
401 NumberFree(LCMb,"execarg"); NumberFree(LCMc,"execarg");
402 goto execargerr;
403 }
404 kGCD = kGCD2;
405 for ( ii = 0; ii < kGCD; ii++ ) GCDbuffer[ii] = GCDbuffer2[ii];
406 }
407 else {
408 kGCD = 1; GCDbuffer[0] = 1;
409 }
410 k = REDLENG(j);
411 if ( k < 0 ) k = -k;
412 r3 += k;
413 while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
414 if ( ( ( LCMbuffer[0] == 1 ) && ( kLCM == 1 ) ) ) {
415 for ( kLCM = 0; kLCM < k; kLCM++ )
416 LCMbuffer[kLCM] = r3[kLCM];
417 }
418 else if ( ( k != 1 ) || ( r3[0] != 1 ) ) {
419 if ( GcdLong(BHEAD LCMbuffer,kLCM,(UWORD *)r3,k,LCMb,&kkLCM) ) {
420 NumberFree(GCDbuffer,"execarg"); NumberFree(GCDbuffer2,"execarg");
421 NumberFree(LCMbuffer,"execarg"); NumberFree(LCMb,"execarg"); NumberFree(LCMc,"execarg");
422 goto execargerr;
423 }
424 DivLong((UWORD *)r3,k,LCMb,kkLCM,LCMb,&kkLCM,LCMc,&jLCM);
425 MulLong(LCMbuffer,kLCM,LCMb,kkLCM,LCMc,&jLCM);
426 for ( kLCM = 0; kLCM < jLCM; kLCM++ )
427 LCMbuffer[kLCM] = LCMc[kLCM];
428 }
429 else {} /* LCM doesn't change */
430 r1 = r2;
431 }
432/*
433 Now put the factor together: GCD/LCM
434*/
435 r3 = (WORD *)(GCDbuffer);
436 if ( kGCD == kLCM ) {
437 for ( jGCD = 0; jGCD < kGCD; jGCD++ )
438 r3[jGCD+kGCD] = LCMbuffer[jGCD];
439 k = kGCD;
440 }
441 else if ( kGCD > kLCM ) {
442 for ( jGCD = 0; jGCD < kLCM; jGCD++ )
443 r3[jGCD+kGCD] = LCMbuffer[jGCD];
444 for ( jGCD = kLCM; jGCD < kGCD; jGCD++ )
445 r3[jGCD+kGCD] = 0;
446 k = kGCD;
447 }
448 else {
449 for ( jGCD = kGCD; jGCD < kLCM; jGCD++ )
450 r3[jGCD] = 0;
451 for ( jGCD = 0; jGCD < kLCM; jGCD++ )
452 r3[jGCD+kLCM] = LCMbuffer[jGCD];
453 k = kLCM;
454 }
455
456 j = 2*k+1;
457/*
458 Now we have to correct the overall factor
459*/
460 if ( scale && ( factor == 0 || *factor > 0 ) )
461 goto ScaledVariety;
462/*
463 The if was added 28-nov-2012 to give MakeInteger also
464 the (0) option.
465*/
466 if ( scale && ( factor == 0 || *factor ) ) {
467 size = term[*term-1];
468 size = REDLENG(size);
469 if ( MulRat(BHEAD (UWORD *)rstop,size,(UWORD *)r3,k,
470 (UWORD *)rstop,&size) ) goto execargerr;
471 size = INCLENG(size);
472 k = size < 0 ? -size: size;
473 rstop[k-1] = size*sign;
474 *term = (WORD)(rstop - term) + k;
475 }
476 }
477 else {
478 if ( factor && *factor >= 1 ) {
479 r4 = r + *r;
480 r1 = r + ARGHEAD;
481 while ( r1 < r4 ) {
482 r2 = r1 + *r1;
483 r3 = r2 - ABS(r2[-1]);
484 j = r3 - r1;
485 r5 = factor;
486 if ( j != *r5 ) { r1 = r2; continue; }
487 r5++; r6 = r1+1;
488 while ( --j > 0 ) {
489 if ( *r5 != *r6 ) break;
490 r5++; r6++;
491 }
492 if ( j > 0 ) { r1 = r2; continue; }
493 break;
494 }
495 if ( r1 >= r4 ) continue;
496 }
497 else {
498 r1 = r + ARGHEAD;
499 r2 = r1 + *r1;
500 r3 = r2 - ABS(r2[-1]);
501 }
502 if ( *r3 == 1 && r3[1] == 1 ) {
503 if ( r2[-1] == 3 ) continue;
504 if ( r2[-1] == -3 && type == TYPENORM3 ) continue;
505 }
506 action = 1;
507 v[2] |= DIRTYFLAG;
508 j = r2[-1];
509 k = REDLENG(j);
510 if ( j < 0 ) j = -j;
511 if ( type == TYPENORM && scale && ( factor == 0 || *factor ) ) {
512/*
513 Now we correct the overall factor
514*/
515ScaledVariety:;
516 size = term[*term-1];
517 size = REDLENG(size);
518 if ( scale > 0 ) {
519 for ( jj = 0; jj < scale; jj++ ) {
520 if ( MulRat(BHEAD (UWORD *)rstop,size,(UWORD *)r3,k,
521 (UWORD *)rstop,&size) ) goto execargerr;
522 }
523 }
524 else {
525 for ( jj = 0; jj > scale; jj-- ) {
526 if ( DivRat(BHEAD (UWORD *)rstop,size,(UWORD *)r3,k,
527 (UWORD *)rstop,&size) ) goto execargerr;
528 }
529 }
530 size = INCLENG(size);
531 k = size < 0 ? -size: size;
532 rstop[k-1] = size*sign;
533 *term = (WORD)(rstop - term) + k;
534 }
535 }
536/*
537 We generate a statement for adapting all terms in the
538 argument successively
539*/
540 r4 = AddRHS(AT.ebufnum,1);
541 while ( (r4+j+12) > CC->Top ) r4 = DoubleCbuffer(AT.ebufnum,r4,3);
542 *r4++ = j+1;
543 i = (j-1)/2; /* was (j-1)*2 ????? 17-oct-2017 */
544 for ( k = 0; k < i; k++ ) *r4++ = r3[i+k];
545 for ( k = 0; k < i; k++ ) *r4++ = r3[k];
546 if ( ( type == TYPENORM3 ) || ( type == TYPENORM4 ) ) *r4++ = j*sign;
547 else *r4++ = r3[j-1];
548 *r4++ = 0;
549 CC->rhs[CC->numrhs+1] = r4;
550 CC->Pointer = r4;
551 AT.mulpat[5] = CC->numrhs;
552 AT.mulpat[7] = AT.ebufnum;
553 }
554 else if ( type == TYPEARGTOEXTRASYMBOL ) {
555 WORD n;
556 if ( r[0] < 0 ) {
557 /* The argument is in the fast notation. */
558 WORD tmp[MaX(9,FUNHEAD+5)];
559 switch ( r[0] ) {
560 case -SNUMBER:
561 if ( r[1] == 0 ) {
562 tmp[0] = 0;
563 }
564 else {
565 tmp[0] = 4;
566 tmp[1] = ABS(r[1]);
567 tmp[2] = 1;
568 tmp[3] = r[1] > 0 ? 3 : -3;
569 tmp[4] = 0;
570 }
571 break;
572 case -SYMBOL:
573 tmp[0] = 8;
574 tmp[1] = SYMBOL;
575 tmp[2] = 4;
576 tmp[3] = r[1];
577 tmp[4] = 1;
578 tmp[5] = 1;
579 tmp[6] = 1;
580 tmp[7] = 3;
581 tmp[8] = 0;
582 break;
583 case -INDEX:
584 case -VECTOR:
585 case -MINVECTOR:
586 tmp[0] = 7;
587 tmp[1] = INDEX;
588 tmp[2] = 3;
589 tmp[3] = r[1];
590 tmp[4] = 1;
591 tmp[5] = 1;
592 tmp[6] = r[0] != -MINVECTOR ? 3 : -3;
593 tmp[7] = 0;
594 break;
595 default:
596 if ( r[0] <= -FUNCTION ) {
597 tmp[0] = FUNHEAD+4;
598 tmp[1] = -r[0];
599 tmp[2] = FUNHEAD;
600 ZeroFillRange(tmp,3,1+FUNHEAD);
601 tmp[FUNHEAD+1] = 1;
602 tmp[FUNHEAD+2] = 1;
603 tmp[FUNHEAD+3] = 3;
604 tmp[FUNHEAD+4] = 0;
605 break;
606 }
607 else {
608 MLOCK(ErrorMessageLock);
609 MesPrint("Unknown fast notation found (TYPEARGTOEXTRASYMBOL)");
610 MUNLOCK(ErrorMessageLock);
611 return(-1);
612 }
613 }
614 n = FindSubexpression(tmp);
615 }
616 else {
617 /*
618 * NOTE: writing to r[r[0]] is legal. As long as we work
619 * in a part of the term, at least the coefficient of
620 * the term must follow.
621 */
622 WORD old_rr0 = r[r[0]];
623 r[r[0]] = 0; /* zero-terminated */
624 n = FindSubexpression(r+ARGHEAD);
625 r[r[0]] = old_rr0;
626 }
627 /* Put the new argument in the work space. */
628 if ( AT.WorkPointer+2 > AT.WorkTop ) {
629 MLOCK(ErrorMessageLock);
630 MesWork();
631 MUNLOCK(ErrorMessageLock);
632 return(-1);
633 }
634 r1 = AT.WorkPointer;
635 if ( scale ) { /* means "tonumber" */
636 r1[0] = -SNUMBER;
637 r1[1] = n;
638 }
639 else {
640 r1[0] = -SYMBOL;
641 r1[1] = MAXVARIABLES-n;
642 }
643 /* We need r2, r3, m and k to shift the data. */
644 r2 = r + (r[0] > 0 ? r[0] : r[0] <= -FUNCTION ? 1 : 2);
645 r3 = r;
646 m = r1+ARGHEAD+2;
647 k = 2;
648 goto do_shift;
649 }
650 r3 = r;
651 AR.DeferFlag = 0;
652 if ( *r > 0 ) {
653 NewSort(BHEAD0);
654 action = 1;
655 r2 = r + *r;
656 r += ARGHEAD;
657 while ( r < r2 ) { /* Sum over the terms */
658 m = AT.WorkPointer;
659 j = *r;
660 while ( --j >= 0 ) *m++ = *r++;
661 r1 = AT.WorkPointer;
662 AT.WorkPointer = m;
663/*
664 What to do with dummy indices?
665*/
666 if ( type == TYPENORM || type == TYPENORM2 || type == TYPENORM3 || type == TYPENORM4 ) {
667 if ( MultDo(BHEAD r1,AT.mulpat) ) goto execargerr;
668 AT.WorkPointer = r1 + *r1;
669 }
670 if ( Generator(BHEAD r1,level) ) goto execargerr;
671 AT.WorkPointer = r1;
672 }
673 }
674 else {
675 r2 = r + (( *r <= -FUNCTION ) ? 1:2);
676 r1 = AT.WorkPointer;
677 ToGeneral(r,r1,0);
678 m = r1 + ARGHEAD;
679 AT.WorkPointer = r1 + *r1;
680 NewSort(BHEAD0);
681 action = 1;
682/*
683 What to do with dummy indices?
684*/
685 if ( type == TYPENORM || type == TYPENORM2 || type == TYPENORM3 || type == TYPENORM4 ) {
686 if ( MultDo(BHEAD m,AT.mulpat) ) goto execargerr;
687 AT.WorkPointer = m + *m;
688 }
689 if ( (*m != 0 ) && Generator(BHEAD m,level) ) goto execargerr;
690 AT.WorkPointer = r1;
691 }
692 if ( EndSort(BHEAD AT.WorkPointer+ARGHEAD,1) < 0 ) goto execargerr;
693 AR.DeferFlag = olddefer;
694/*
695 Now shift the sorted entity over the old argument.
696*/
697 m = AT.WorkPointer+ARGHEAD;
698 while ( *m ) m += *m;
699 k = WORDDIF(m,AT.WorkPointer);
700 *AT.WorkPointer = k;
701 AT.WorkPointer[1] = 0;
702 if ( ToFast(AT.WorkPointer,AT.WorkPointer) ) {
703 if ( *AT.WorkPointer <= -FUNCTION ) k = 1;
704 else k = 2;
705 }
706do_shift:
707 if ( *r3 > 0 ) j = k - *r3;
708 else if ( *r3 <= -FUNCTION ) j = k - 1;
709 else j = k - 2;
710
711 t[1] += j;
712 action = 1;
713 v[2] |= DIRTYFLAG;
714 if ( j > 0 ) {
715 r = m + j;
716 while ( m > AT.WorkPointer ) *--r = *--m;
717 AT.WorkPointer = r;
718 m = term + *term;
719 r = m + j;
720 while ( m > r2 ) *--r = *--m;
721 }
722 else if ( j < 0 ) {
723 r = r2 + j;
724 r1 = term + *term;
725 while ( r2 < r1 ) *r++ = *r2++;
726 }
727 r = r3;
728 m = AT.WorkPointer;
729 NCOPY(r,m,k);
730 *term += j;
731 rstop += j;
732 CC->numrhs = oldnumrhs;
733 CC->Pointer = CC->Buffer + oldcpointer;
734 }
735 }
736 t += t[1];
737 }
738/*
739 If TYPENORM4, we allocated Number buffers before the above while loop. Free them.
740*/
741 if ( type == TYPENORM4 ) {
742 NumberFree(GCDbuffer,"execarg");
743 NumberFree(GCDbuffer2,"execarg");
744 NumberFree(LCMbuffer,"execarg");
745 NumberFree(LCMb,"execarg"); NumberFree(LCMc,"execarg");
746 }
747 if ( didpolyratfun ) {
748 PolyFunDirty(BHEAD term);
749 didpolyratfun = 0;
750 }
751/*
752 #] Argument detection :
753 #[ SplitArg : + varieties
754*/
755 if ( ( type == TYPESPLITARG || type == TYPESPLITARG2
756 || type == TYPESPLITFIRSTARG || type == TYPESPLITLASTARG ) &&
757 AT.pWorkPointer > oldppointer ) {
758 t = term+1;
759 r1 = AT.WorkPointer + 1;
760 lp = oldppointer;
761 while ( t < rstop ) {
762 if ( lp < AT.pWorkPointer && t == AT.pWorkSpace[lp] ) {
763 v = t;
764 m = t + FUNHEAD;
765 r = t + t[1];
766 r2 = r1; while ( t < m ) *r1++ = *t++;
767 while ( m < r ) {
768 t = m;
769 NEXTARG(m)
770 if ( lp >= AT.pWorkPointer || t != AT.pWorkSpace[lp+1] ) {
771 if ( *t > 0 ) t[1] = 0;
772 while ( t < m ) *r1++ = *t++;
773 continue;
774 }
775/*
776 Now we have a nontrivial argument that should be done.
777*/
778 lp += 2;
779 action = 1;
780 v[2] |= DIRTYFLAG;
781 r3 = t + *t;
782 t += ARGHEAD;
783 if ( type == TYPESPLITFIRSTARG ) {
784 r4 = r1; r5 = t; r7 = oldwork;
785 *r1++ = *t + ARGHEAD;
786 for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
787 j = 0;
788 while ( t < r3 ) {
789 i = *t;
790 if ( j == 0 ) {
791 NCOPY(r7,t,i)
792 j++;
793 }
794 else {
795 NCOPY(r1,t,i)
796 }
797 }
798 *r4 = r1 - r4;
799 if ( j ) {
800 if ( ToFast(r4,r4) ) {
801 r1 = r4;
802 if ( *r1 > -FUNCTION ) r1++;
803 r1++;
804 }
805 r7 = oldwork;
806 while ( --j >= 0 ) {
807 r4 = r1; i = *r7;
808 *r1++ = i+ARGHEAD; *r1++ = 0;
809 FILLARG(r1);
810 NCOPY(r1,r7,i)
811 if ( ToFast(r4,r4) ) {
812 r1 = r4;
813 if ( *r1 > -FUNCTION ) r1++;
814 r1++;
815 }
816 }
817 }
818 t = r3;
819 }
820 else if ( type == TYPESPLITLASTARG ) {
821 r4 = r1; r5 = t; r7 = oldwork;
822 *r1++ = *t + ARGHEAD;
823 for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
824 j = 0;
825 while ( t < r3 ) {
826 i = *t;
827 if ( t+i >= r3 ) {
828 NCOPY(r7,t,i)
829 j++;
830 }
831 else {
832 NCOPY(r1,t,i)
833 }
834 }
835 *r4 = r1 - r4;
836 if ( j ) {
837 if ( ToFast(r4,r4) ) {
838 r1 = r4;
839 if ( *r1 > -FUNCTION ) r1++;
840 r1++;
841 }
842 r7 = oldwork;
843 while ( --j >= 0 ) {
844 r4 = r1; i = *r7;
845 *r1++ = i+ARGHEAD; *r1++ = 0;
846 FILLARG(r1);
847 NCOPY(r1,r7,i)
848 if ( ToFast(r4,r4) ) {
849 r1 = r4;
850 if ( *r1 > -FUNCTION ) r1++;
851 r1++;
852 }
853 }
854 }
855 t = r3;
856 }
857 else if ( factor == 0 || ( type == TYPESPLITARG2 && *factor == 0 ) ) {
858 while ( t < r3 ) {
859 r4 = r1;
860 *r1++ = *t + ARGHEAD;
861 for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
862 i = *t;
863 while ( --i >= 0 ) *r1++ = *t++;
864 if ( ToFast(r4,r4) ) {
865 r1 = r4;
866 if ( *r1 > -FUNCTION ) r1++;
867 r1++;
868 }
869 }
870 }
871 else if ( type == TYPESPLITARG2 ) {
872/*
873 Here we better put the pattern matcher at work?
874 Remember: there are no wildcards.
875*/
876 WORD *oRepFunList = AN.RepFunList;
877 WORD *oWildMask = AT.WildMask, *oWildValue = AN.WildValue;
878 AN.WildValue = AT.locwildvalue; AT.WildMask = AT.locwildvalue+2;
879 AN.NumWild = 0;
880 r4 = r1; r5 = t; r7 = oldwork;
881 *r1++ = *t + ARGHEAD;
882 for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
883 j = 0;
884 while ( t < r3 ) {
885 AN.UseFindOnly = 0; oldwork2 = AT.WorkPointer;
886 AN.RepFunList = r1;
887 AT.WorkPointer = r1+AN.RepFunNum+2;
888 i = *t;
889 if ( FindRest(BHEAD t,factor) &&
890 ( AN.UsedOtherFind || FindOnce(BHEAD t,factor) ) ) {
891 NCOPY(r7,t,i)
892 j++;
893 }
894 else if ( factor[0] == FUNHEAD+1 && factor[1] >= FUNCTION ) {
895 WORD *rr1 = t+1, *rr2 = t+i;
896 rr2 -= ABS(rr2[-1]);
897 while ( rr1 < rr2 ) {
898 if ( *rr1 == factor[1] ) break;
899 rr1 += rr1[1];
900 }
901 if ( rr1 < rr2 ) {
902 NCOPY(r7,t,i)
903 j++;
904 }
905 else {
906 NCOPY(r1,t,i)
907 }
908 }
909 else {
910 NCOPY(r1,t,i)
911 }
912 AT.WorkPointer = oldwork2;
913 }
914 AN.RepFunList = oRepFunList;
915 *r4 = r1 - r4;
916 if ( j ) {
917 if ( ToFast(r4,r4) ) {
918 r1 = r4;
919 if ( *r1 > -FUNCTION ) r1++;
920 r1++;
921 }
922 r7 = oldwork;
923 while ( --j >= 0 ) {
924 r4 = r1; i = *r7;
925 *r1++ = i+ARGHEAD; *r1++ = 0;
926 FILLARG(r1);
927 NCOPY(r1,r7,i)
928 if ( ToFast(r4,r4) ) {
929 r1 = r4;
930 if ( *r1 > -FUNCTION ) r1++;
931 r1++;
932 }
933 }
934 }
935 t = r3;
936 AT.WildMask = oWildMask; AN.WildValue = oWildValue;
937 }
938 else {
939/*
940 This code deals with splitting off a single term
941*/
942 r4 = r1; r5 = t;
943 *r1++ = *t + ARGHEAD;
944 for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
945 j = 0;
946 while ( t < r3 ) {
947 r6 = t + *t; r6 -= ABS(r6[-1]);
948 if ( (r6 - t) == *factor ) {
949 k = *factor - 1;
950 for ( ; k > 0; k-- ) {
951 if ( t[k] != factor[k] ) break;
952 }
953 if ( k <= 0 ) {
954 j = r3 - t; t += *t; continue;
955 }
956 }
957 else if ( (r6 - t) == 1 && *factor == 0 ) {
958 j = r3 - t; t += *t; continue;
959 }
960 i = *t;
961 NCOPY(r1,t,i)
962 }
963 *r4 = r1 - r4;
964 if ( j ) {
965 if ( ToFast(r4,r4) ) {
966 r1 = r4;
967 if ( *r1 > -FUNCTION ) r1++;
968 r1++;
969 }
970 t = r3 - j;
971 r4 = r1;
972 *r1++ = *t + ARGHEAD;
973 for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
974 i = *t;
975 while ( --i >= 0 ) *r1++ = *t++;
976 if ( ToFast(r4,r4) ) {
977 r1 = r4;
978 if ( *r1 > -FUNCTION ) r1++;
979 r1++;
980 }
981 }
982 t = r3;
983 }
984 }
985 r2[1] = r1 - r2;
986 }
987 else {
988 r = t + t[1];
989 while ( t < r ) *r1++ = *t++;
990 }
991 }
992 r = term + *term;
993 while ( t < r ) *r1++ = *t++;
994 m = AT.WorkPointer;
995 i = m[0] = r1 - m;
996 t = term;
997 while ( --i >= 0 ) *t++ = *m++;
998 if ( AT.WorkPointer < m ) AT.WorkPointer = m;
999 }
1000/*
1001 #] SplitArg :
1002 #[ FACTARG :
1003*/
1004 if ( ( type == TYPEFACTARG || type == TYPEFACTARG2 ) &&
1005 AT.pWorkPointer > oldppointer ) {
1006 t = term+1;
1007 r1 = AT.WorkPointer + 1;
1008 lp = oldppointer;
1009 while ( t < rstop ) {
1010 if ( lp < AT.pWorkPointer && AT.pWorkSpace[lp] == t ) {
1011 v = t;
1012 m = t + FUNHEAD;
1013 r = t + t[1];
1014 r2 = r1; while ( t < m ) *r1++ = *t++;
1015 while ( m < r ) {
1016 rr = t = m;
1017 NEXTARG(m)
1018 if ( lp >= AT.pWorkPointer || AT.pWorkSpace[lp+1] != t ) {
1019 if ( *t > 0 ) t[1] = 0;
1020 while ( t < m ) *r1++ = *t++;
1021 continue;
1022 }
1023/*
1024 Now we have a nontrivial argument that should be studied.
1025 Try to find common factors.
1026*/
1027 lp += 2;
1028 if ( *t < 0 ) {
1029 if ( factor && ( *factor == 0 && *t == -SNUMBER ) ) {
1030 *r1++ = *t++;
1031 if ( *t == 0 ) *r1++ = *t++;
1032 else { *r1++ = 1; t++; }
1033 continue;
1034 }
1035 else if ( factor && *factor == 4 && factor[2] == 1 ) {
1036 if ( *t == -SNUMBER ) {
1037 if ( factor[3] < 0 || t[1] >= 0 ) {
1038 while ( t < m ) *r1++ = *t++;
1039 }
1040 else {
1041 *r1++ = -SNUMBER; *r1++ = -1;
1042 *r1++ = *t++; *r1++ = -*t++;
1043 }
1044 }
1045 else {
1046 while ( t < m ) *r1++ = *t++;
1047 *r1++ = -SNUMBER; *r1++ = 1;
1048 }
1049 continue;
1050 }
1051 else if ( *t == -MINVECTOR ) {
1052 if ( AC.OldFactArgFlag == NEWFACTARG ) {
1053 *r1++ = -SNUMBER; *r1++ = -1;
1054 *r1++ = -VECTOR; t++; *r1++ = *t++;
1055 }
1056 else {
1057 *r1++ = -VECTOR; t++; *r1++ = *t++;
1058 *r1++ = -SNUMBER; *r1++ = -1;
1059 *r1++ = -SNUMBER; *r1++ = 1;
1060 }
1061 continue;
1062 }
1063 }
1064/*
1065 Now we have a nontrivial argument
1066*/
1067 r3 = t + *t;
1068 t += ARGHEAD; r5 = t; /* Store starting point */
1069 /* We have terms from r5 to r3 */
1070 if ( r5+*r5 == r3 && factor ) { /* One term only */
1071 if ( *factor == 0 ) {
1072 GETSTOP(t,r6);
1073 r9 = r1; *r1++ = 0; *r1++ = 1;
1074 FILLARG(r1);
1075 *r1++ = (r6-t)+3; t++;
1076 while ( t < r6 ) *r1++ = *t++;
1077 *r1++ = 1; *r1++ = 1; *r1++ = 3;
1078 *r9 = r1-r9;
1079 if ( ToFast(r9,r9) ) {
1080 if ( *r9 <= -FUNCTION ) r1 = r9+1;
1081 else r1 = r9+2;
1082 }
1083 t = r3; continue;
1084 }
1085 if ( factor[0] == 4 && factor[2] == 1 ) {
1086 GETSTOP(t,r6);
1087 r7 = r1; *r1++ = (r6-t)+3+ARGHEAD; *r1++ = 0;
1088 FILLARG(r1);
1089 *r1++ = (r6-t)+3; t++;
1090 while ( t < r6 ) *r1++ = *t++;
1091 *r1++ = 1; *r1++ = 1; *r1++ = 3;
1092 if ( ToFast(r7,r7) ) {
1093 if ( *r7 <= -FUNCTION ) r1 = r7+1;
1094 else r1 = r7+2;
1095 }
1096 if ( r3[-1] < 0 && factor[3] > 0 ) {
1097 *r1++ = -SNUMBER; *r1++ = -1;
1098 if ( r3[-1] == -3 && r3[-2] == 1
1099 && ( r3[-3] & MAXPOSITIVE ) == r3[-3] ) {
1100 *r1++ = -SNUMBER; *r1++ = r3[-3];
1101 }
1102 else {
1103 *r1++ = (r3-r6)+1+ARGHEAD;
1104 *r1++ = 0;
1105 FILLARG(r1);
1106 *r1++ = (r3-r6+1);
1107 while ( t < r3 ) *r1++ = *t++;
1108 r1[-1] = -r1[-1];
1109 }
1110 }
1111 else {
1112 if ( ( r3[-1] == -3 || r3[-1] == 3 )
1113 && r3[-2] == 1
1114 && ( r3[-3] & MAXPOSITIVE ) == r3[-3] ) {
1115 *r1++ = -SNUMBER; *r1++ = r3[-3];
1116 if ( r3[-1] < 0 ) r1[-1] = - r1[-1];
1117 }
1118 else {
1119 *r1++ = (r3-r6)+1+ARGHEAD;
1120 *r1++ = 0;
1121 FILLARG(r1);
1122 *r1++ = (r3-r6+1);
1123 while ( t < r3 ) *r1++ = *t++;
1124 }
1125 }
1126 t = r3; continue;
1127 }
1128 }
1129/*
1130 Now we take the first term and look for its pieces
1131 inside the other terms.
1132
1133 It is at this point that a more general factorization
1134 routine could take over (allowing for writing the output
1135 properly of course).
1136*/
1137 if ( AC.OldFactArgFlag == NEWFACTARG ) {
1138 if ( factor == 0 ) {
1139 WORD *oldworkpointer2 = AT.WorkPointer;
1140 AT.WorkPointer = r1 + AM.MaxTer+FUNHEAD;
1141 if ( ArgFactorize(BHEAD t-ARGHEAD,r1) < 0 ) {
1142 MesCall("ExecArg");
1143 return(-1);
1144 }
1145 AT.WorkPointer = oldworkpointer2;
1146 t = r3;
1147 while ( *r1 ) { NEXTARG(r1) }
1148 }
1149 else {
1150 rnext = t + *t;
1151 GETSTOP(t,r6);
1152 t++;
1153 t = r5; pow = 1;
1154 while ( t < r3 ) {
1155 t += *t; if ( t[-1] > 0 ) { pow = 0; break; }
1156 }
1157/*
1158 We have to add here the code for computing the GCD
1159 and to divide it out.
1160
1161 #[ Numerical factor :
1162*/
1163 t = r5;
1164 EAscrat = (UWORD *)(TermMalloc("execarg"));
1165 if ( t + *t == r3 ) {
1166 if ( factor == 0 || *factor > 2 ) {
1167 if ( pow > 0 ) {
1168 *r1++ = -SNUMBER; *r1++ = -1;
1169 t = r5;
1170 while ( t < r3 ) {
1171 t += *t; t[-1] = -t[-1];
1172 }
1173 }
1174 t = rr; *r1++ = *t++; *r1++ = 1; t++;
1175 COPYARG(r1,t);
1176 while ( t < m ) *r1++ = *t++;
1177 }
1178 }
1179 else {
1180 GETSTOP(t,r6);
1181 ngcd = t[t[0]-1];
1182 i = abs(ngcd)-1;
1183 while ( --i >= 0 ) EAscrat[i] = r6[i];
1184 t += *t;
1185 while ( t < r3 ) {
1186 GETSTOP(t,r6);
1187 i = t[t[0]-1];
1188 if ( AccumGCD(BHEAD EAscrat,&ngcd,(UWORD *)r6,i) ) goto execargerr;
1189 if ( ngcd == 3 && EAscrat[0] == 1 && EAscrat[1] == 1 ) break;
1190 t += *t;
1191 }
1192/*
1193 if ( ngcd != 3 || EAscrat[0] != 1 || EAscrat[1] != 1 )
1194*/
1195 {
1196 if ( pow ) ngcd = -ngcd;
1197 t = r5; r9 = r1; *r1++ = t[-ARGHEAD]; *r1++ = 1;
1198 FILLARG(r1); ngcd = REDLENG(ngcd);
1199 while ( t < r3 ) {
1200 GETSTOP(t,r6);
1201 r7 = t; r8 = r1;
1202 while ( r7 < r6) *r1++ = *r7++;
1203 t += *t;
1204 i = REDLENG(t[-1]);
1205 if ( DivRat(BHEAD (UWORD *)r6,i,EAscrat,ngcd,(UWORD *)r1,&nq) ) goto execargerr;
1206 nq = INCLENG(nq);
1207 i = ABS(nq)-1;
1208 r1 += i; *r1++ = nq; *r8 = r1-r8;
1209 }
1210 *r9 = r1-r9;
1211 ngcd = INCLENG(ngcd);
1212 i = ABS(ngcd)-1;
1213 if ( factor && *factor == 0 ) {}
1214 else if ( ( factor && factor[0] == 4 && factor[2] == 1
1215 && factor[3] == -3 ) || pow == 0 ) {
1216 r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1217 FILLARG(r1); *r1++ = i+2;
1218 for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1219 *r1++ = ngcd;
1220 if ( ToFast(r9,r9) ) r1 = r9+2;
1221 }
1222 else if ( factor && factor[0] == 4 && factor[2] == 1
1223 && factor[3] > 0 && pow ) {
1224 if ( ngcd < 0 ) ngcd = -ngcd;
1225 *r1++ = -SNUMBER; *r1++ = -1;
1226 r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1227 FILLARG(r1); *r1++ = i+2;
1228 for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1229 *r1++ = ngcd;
1230 if ( ToFast(r9,r9) ) r1 = r9+2;
1231 }
1232 else {
1233 if ( ngcd < 0 ) ngcd = -ngcd;
1234 if ( pow ) { *r1++ = -SNUMBER; *r1++ = -1; }
1235 if ( ngcd != 3 || EAscrat[0] != 1 || EAscrat[1] != 1 ) {
1236 r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1237 FILLARG(r1); *r1++ = i+2;
1238 for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1239 *r1++ = ngcd;
1240 if ( ToFast(r9,r9) ) r1 = r9+2;
1241 }
1242 }
1243 }
1244/*
1245 #] Numerical factor :
1246 else {
1247onetermnew:;
1248
1249 if ( factor == 0 || *factor > 2 ) {
1250 if ( pow > 0 ) {
1251 *r1++ = -SNUMBER; *r1++ = -1;
1252 t = r5;
1253 while ( t < r3 ) {
1254 t += *t; t[-1] = -t[-1];
1255 }
1256 }
1257 t = rr; *r1++ = *t++; *r1++ = 1; t++;
1258 COPYARG(r1,t);
1259 while ( t < m ) *r1++ = *t++;
1260 }
1261 }
1262onetermnew:;
1263*/
1264 }
1265 TermFree(EAscrat,"execarg");
1266 }
1267 }
1268 else { /* AC.OldFactArgFlag is ON */
1269 {
1270 WORD *mnext, ncom;
1271 rnext = t + *t;
1272 GETSTOP(t,r6);
1273 t++;
1274 if ( factor == 0 ) {
1275 while ( t < r6 ) {
1276/*
1277 #[ SYMBOL :
1278*/
1279 if ( *t == SYMBOL ) {
1280 r7 = t; r8 = t + t[1]; t += 2;
1281 while ( t < r8 ) {
1282 pow = t[1];
1283 mm = rnext;
1284 while ( mm < r3 ) {
1285 mnext = mm + *mm;
1286 GETSTOP(mm,mstop); mm++;
1287 while ( mm < mstop ) {
1288 if ( *mm != SYMBOL ) mm += mm[1];
1289 else break;
1290 }
1291 if ( *mm == SYMBOL ) {
1292 mstop = mm + mm[1]; mm += 2;
1293 while ( *mm != *t && mm < mstop ) mm += 2;
1294 if ( mm >= mstop ) pow = 0;
1295 else if ( pow > 0 && mm[1] > 0 ) {
1296 if ( mm[1] < pow ) pow = mm[1];
1297 }
1298 else if ( pow < 0 && mm[1] < 0 ) {
1299 if ( mm[1] > pow ) pow = mm[1];
1300 }
1301 else pow = 0;
1302 }
1303 else pow = 0;
1304 if ( pow == 0 ) break;
1305 mm = mnext;
1306 }
1307 if ( pow == 0 ) { t += 2; continue; }
1308/*
1309 We have a factor
1310*/
1311 action = 1; i = pow;
1312 if ( i > 0 ) {
1313 while ( --i >= 0 ) {
1314 *r1++ = -SYMBOL;
1315 *r1++ = *t;
1316 }
1317 }
1318 else {
1319 while ( i++ < 0 ) {
1320 *r1++ = 8 + ARGHEAD;
1321 for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1322 *r1++ = 8; *r1++ = SYMBOL;
1323 *r1++ = 4; *r1++ = *t; *r1++ = -1;
1324 *r1++ = 1; *r1++ = 1; *r1++ = 3;
1325 }
1326 }
1327/*
1328 Now we have to remove the symbols
1329*/
1330 t[1] -= pow;
1331 mm = rnext;
1332 while ( mm < r3 ) {
1333 mnext = mm + *mm;
1334 GETSTOP(mm,mstop); mm++;
1335 while ( mm < mstop ) {
1336 if ( *mm != SYMBOL ) mm += mm[1];
1337 else break;
1338 }
1339 mstop = mm + mm[1]; mm += 2;
1340 while ( mm < mstop && *mm != *t ) mm += 2;
1341 mm[1] -= pow;
1342 mm = mnext;
1343 }
1344 t += 2;
1345 }
1346 }
1347/*
1348 #] SYMBOL :
1349 #[ DOTPRODUCT :
1350*/
1351 else if ( *t == DOTPRODUCT ) {
1352 r7 = t; r8 = t + t[1]; t += 2;
1353 while ( t < r8 ) {
1354 pow = t[2];
1355 mm = rnext;
1356 while ( mm < r3 ) {
1357 mnext = mm + *mm;
1358 GETSTOP(mm,mstop); mm++;
1359 while ( mm < mstop ) {
1360 if ( *mm != DOTPRODUCT ) mm += mm[1];
1361 else break;
1362 }
1363 if ( *mm == DOTPRODUCT ) {
1364 mstop = mm + mm[1]; mm += 2;
1365 while ( ( *mm != *t || mm[1] != t[1] )
1366 && mm < mstop ) mm += 3;
1367 if ( mm >= mstop ) pow = 0;
1368 else if ( pow > 0 && mm[2] > 0 ) {
1369 if ( mm[2] < pow ) pow = mm[2];
1370 }
1371 else if ( pow < 0 && mm[2] < 0 ) {
1372 if ( mm[2] > pow ) pow = mm[2];
1373 }
1374 else pow = 0;
1375 }
1376 else pow = 0;
1377 if ( pow == 0 ) break;
1378 mm = mnext;
1379 }
1380 if ( pow == 0 ) { t += 3; continue; }
1381/*
1382 We have a factor
1383*/
1384 action = 1; i = pow;
1385 if ( i > 0 ) {
1386 while ( --i >= 0 ) {
1387 *r1++ = 9 + ARGHEAD;
1388 for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1389 *r1++ = 9; *r1++ = DOTPRODUCT;
1390 *r1++ = 5; *r1++ = *t; *r1++ = t[1]; *r1++ = 1;
1391 *r1++ = 1; *r1++ = 1; *r1++ = 3;
1392 }
1393 }
1394 else {
1395 while ( i++ < 0 ) {
1396 *r1++ = 9 + ARGHEAD;
1397 for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1398 *r1++ = 9; *r1++ = DOTPRODUCT;
1399 *r1++ = 5; *r1++ = *t; *r1++ = t[1]; *r1++ = -1;
1400 *r1++ = 1; *r1++ = 1; *r1++ = 3;
1401 }
1402 }
1403/*
1404 Now we have to remove the dotproducts
1405*/
1406 t[2] -= pow;
1407 mm = rnext;
1408 while ( mm < r3 ) {
1409 mnext = mm + *mm;
1410 GETSTOP(mm,mstop); mm++;
1411 while ( mm < mstop ) {
1412 if ( *mm != DOTPRODUCT ) mm += mm[1];
1413 else break;
1414 }
1415 mstop = mm + mm[1]; mm += 2;
1416 while ( mm < mstop && ( *mm != *t
1417 || mm[1] != t[1] ) ) mm += 3;
1418 mm[2] -= pow;
1419 mm = mnext;
1420 }
1421 t += 3;
1422 }
1423 }
1424/*
1425 #] DOTPRODUCT :
1426 #[ DELTA/VECTOR :
1427*/
1428 else if ( *t == DELTA || *t == VECTOR ) {
1429 r7 = t; r8 = t + t[1]; t += 2;
1430 while ( t < r8 ) {
1431 mm = rnext;
1432 pow = 1;
1433 while ( mm < r3 ) {
1434 mnext = mm + *mm;
1435 GETSTOP(mm,mstop); mm++;
1436 while ( mm < mstop ) {
1437 if ( *mm != *r7 ) mm += mm[1];
1438 else break;
1439 }
1440 if ( *mm == *r7 ) {
1441 mstop = mm + mm[1]; mm += 2;
1442 while ( ( *mm != *t || mm[1] != t[1] )
1443 && mm < mstop ) mm += 2;
1444 if ( mm >= mstop ) pow = 0;
1445 }
1446 else pow = 0;
1447 if ( pow == 0 ) break;
1448 mm = mnext;
1449 }
1450 if ( pow == 0 ) { t += 2; continue; }
1451/*
1452 We have a factor
1453*/
1454 action = 1;
1455 *r1++ = 8 + ARGHEAD;
1456 for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1457 *r1++ = 8; *r1++ = *r7;
1458 *r1++ = 4; *r1++ = *t; *r1++ = t[1];
1459 *r1++ = 1; *r1++ = 1; *r1++ = 3;
1460/*
1461 Now we have to remove the delta's/vectors
1462*/
1463 mm = rnext;
1464 while ( mm < r3 ) {
1465 mnext = mm + *mm;
1466 GETSTOP(mm,mstop); mm++;
1467 while ( mm < mstop ) {
1468 if ( *mm != *r7 ) mm += mm[1];
1469 else break;
1470 }
1471 mstop = mm + mm[1]; mm += 2;
1472 while ( mm < mstop && (
1473 *mm != *t || mm[1] != t[1] ) ) mm += 2;
1474 *mm = mm[1] = NOINDEX;
1475 mm = mnext;
1476 }
1477 *t = t[1] = NOINDEX;
1478 t += 2;
1479 }
1480 }
1481/*
1482 #] DELTA/VECTOR :
1483 #[ INDEX :
1484*/
1485 else if ( *t == INDEX ) {
1486 r7 = t; r8 = t + t[1]; t += 2;
1487 while ( t < r8 ) {
1488 mm = rnext;
1489 pow = 1;
1490 while ( mm < r3 ) {
1491 mnext = mm + *mm;
1492 GETSTOP(mm,mstop); mm++;
1493 while ( mm < mstop ) {
1494 if ( *mm != *r7 ) mm += mm[1];
1495 else break;
1496 }
1497 if ( *mm == *r7 ) {
1498 mstop = mm + mm[1]; mm += 2;
1499 while ( *mm != *t
1500 && mm < mstop ) mm++;
1501 if ( mm >= mstop ) pow = 0;
1502 }
1503 else pow = 0;
1504 if ( pow == 0 ) break;
1505 mm = mnext;
1506 }
1507 if ( pow == 0 ) { t++; continue; }
1508/*
1509 We have a factor
1510*/
1511 action = 1;
1512/*
1513 The next looks like an error.
1514 We should have here a VECTOR or INDEX like object
1515
1516 *r1++ = 7 + ARGHEAD;
1517 for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
1518 *r1++ = 7; *r1++ = *r7;
1519 *r1++ = 3; *r1++ = *t;
1520 *r1++ = 1; *r1++ = 1; *r1++ = 3;
1521
1522 Replace this by: (11-apr-2007)
1523*/
1524 if ( *t < 0 ) { *r1++ = -VECTOR; }
1525 else { *r1++ = -INDEX; }
1526 *r1++ = *t;
1527/*
1528 Now we have to remove the index
1529*/
1530 *t = NOINDEX;
1531 mm = rnext;
1532 while ( mm < r3 ) {
1533 mnext = mm + *mm;
1534 GETSTOP(mm,mstop); mm++;
1535 while ( mm < mstop ) {
1536 if ( *mm != *r7 ) mm += mm[1];
1537 else break;
1538 }
1539 mstop = mm + mm[1]; mm += 2;
1540 while ( mm < mstop &&
1541 *mm != *t ) mm += 1;
1542 *mm = NOINDEX;
1543 mm = mnext;
1544 }
1545 t += 1;
1546 }
1547 }
1548/*
1549 #] INDEX :
1550 #[ FUNCTION :
1551*/
1552 else if ( *t >= FUNCTION ) {
1553/*
1554 In the next code we should actually look inside
1555 the DENOMINATOR or EXPONENT for noncommuting objects
1556*/
1557 if ( *t >= FUNCTION &&
1558 functions[*t-FUNCTION].commute == 0 ) ncom = 0;
1559 else ncom = 1;
1560 if ( ncom ) {
1561 mm = r5 + 1;
1562 while ( mm < t && ( *mm == DUMMYFUN
1563 || *mm == DUMMYTEN ) ) mm += mm[1];
1564 if ( mm < t ) { t += t[1]; continue; }
1565 }
1566 mm = rnext; pow = 1;
1567 while ( mm < r3 ) {
1568 mnext = mm + *mm;
1569 GETSTOP(mm,mstop); mm++;
1570 while ( mm < mstop ) {
1571 if ( *mm == *t && mm[1] == t[1] ) {
1572 for ( i = 2; i < t[1]; i++ ) {
1573 if ( mm[i] != t[i] ) break;
1574 }
1575 if ( i >= t[1] )
1576 { mm += mm[1]; goto nextmterm; }
1577 }
1578 if ( ncom && *mm != DUMMYFUN && *mm != DUMMYTEN )
1579 { pow = 0; break; }
1580 mm += mm[1];
1581 }
1582 if ( mm >= mstop ) pow = 0;
1583 if ( pow == 0 ) break;
1584nextmterm: mm = mnext;
1585 }
1586 if ( pow == 0 ) { t += t[1]; continue; }
1587/*
1588 Copy the function
1589*/
1590 action = 1;
1591 *r1++ = t[1] + 4 + ARGHEAD;
1592 for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
1593 *r1++ = t[1] + 4;
1594 for ( i = 0; i < t[1]; i++ ) *r1++ = t[i];
1595 *r1++ = 1; *r1++ = 1; *r1++ = 3;
1596/*
1597 Now we have to take out the functions
1598*/
1599 mm = rnext;
1600 while ( mm < r3 ) {
1601 mnext = mm + *mm;
1602 GETSTOP(mm,mstop); mm++;
1603 while ( mm < mstop ) {
1604 if ( *mm == *t && mm[1] == t[1] ) {
1605 for ( i = 2; i < t[1]; i++ ) {
1606 if ( mm[i] != t[i] ) break;
1607 }
1608 if ( i >= t[1] ) {
1609 if ( functions[*t-FUNCTION].spec > 0 )
1610 *mm = DUMMYTEN;
1611 else
1612 *mm = DUMMYFUN;
1613 mm += mm[1];
1614 goto nextterm;
1615 }
1616 }
1617 mm += mm[1];
1618 }
1619nextterm: mm = mnext;
1620 }
1621 if ( functions[*t-FUNCTION].spec > 0 )
1622 *t = DUMMYTEN;
1623 else
1624 *t = DUMMYFUN;
1625 action = 1;
1626 v[2] = DIRTYFLAG;
1627 t += t[1];
1628 }
1629/*
1630 #] FUNCTION :
1631*/
1632 else {
1633 t += t[1];
1634 }
1635 }
1636 }
1637 t = r5; pow = 1;
1638 while ( t < r3 ) {
1639 t += *t; if ( t[-1] > 0 ) { pow = 0; break; }
1640 }
1641/*
1642 We have to add here the code for computing the GCD
1643 and to divide it out.
1644*/
1645/*
1646 #[ Numerical factor :
1647*/
1648 t = r5;
1649 EAscrat = (UWORD *)(TermMalloc("execarg"));
1650 if ( t + *t == r3 ) goto oneterm;
1651 GETSTOP(t,r6);
1652 ngcd = t[t[0]-1];
1653 i = abs(ngcd)-1;
1654 while ( --i >= 0 ) EAscrat[i] = r6[i];
1655 t += *t;
1656 while ( t < r3 ) {
1657 GETSTOP(t,r6);
1658 i = t[t[0]-1];
1659 if ( AccumGCD(BHEAD EAscrat,&ngcd,(UWORD *)r6,i) ) goto execargerr;
1660 if ( ngcd == 3 && EAscrat[0] == 1 && EAscrat[1] == 1 ) break;
1661 t += *t;
1662 }
1663 if ( ngcd != 3 || EAscrat[0] != 1 || EAscrat[1] != 1 ) {
1664 if ( pow ) ngcd = -ngcd;
1665 t = r5; r9 = r1; *r1++ = t[-ARGHEAD]; *r1++ = 1;
1666 FILLARG(r1); ngcd = REDLENG(ngcd);
1667 while ( t < r3 ) {
1668 GETSTOP(t,r6);
1669 r7 = t; r8 = r1;
1670 while ( r7 < r6) *r1++ = *r7++;
1671 t += *t;
1672 i = REDLENG(t[-1]);
1673 if ( DivRat(BHEAD (UWORD *)r6,i,EAscrat,ngcd,(UWORD *)r1,&nq) ) goto execargerr;
1674 nq = INCLENG(nq);
1675 i = ABS(nq)-1;
1676 r1 += i; *r1++ = nq; *r8 = r1-r8;
1677 }
1678 *r9 = r1-r9;
1679 ngcd = INCLENG(ngcd);
1680 i = ABS(ngcd)-1;
1681 if ( factor && *factor == 0 ) {}
1682 else if ( ( factor && factor[0] == 4 && factor[2] == 1
1683 && factor[3] == -3 ) || pow == 0 ) {
1684 r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1685 FILLARG(r1); *r1++ = i+2;
1686 for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1687 *r1++ = ngcd;
1688 if ( ToFast(r9,r9) ) r1 = r9+2;
1689 }
1690 else if ( factor && factor[0] == 4 && factor[2] == 1
1691 && factor[3] > 0 && pow ) {
1692 if ( ngcd < 0 ) ngcd = -ngcd;
1693 *r1++ = -SNUMBER; *r1++ = -1;
1694 r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1695 FILLARG(r1); *r1++ = i+2;
1696 for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1697 *r1++ = ngcd;
1698 if ( ToFast(r9,r9) ) r1 = r9+2;
1699 }
1700 else {
1701 if ( ngcd < 0 ) ngcd = -ngcd;
1702 if ( pow ) { *r1++ = -SNUMBER; *r1++ = -1; }
1703 if ( ngcd != 3 || EAscrat[0] != 1 || EAscrat[1] != 1 ) {
1704 r9 = r1; *r1++ = ARGHEAD+2+i; *r1++ = 0;
1705 FILLARG(r1); *r1++ = i+2;
1706 for ( j = 0; j < i; j++ ) *r1++ = EAscrat[j];
1707 *r1++ = ngcd;
1708 if ( ToFast(r9,r9) ) r1 = r9+2;
1709 }
1710 }
1711 }
1712/*
1713 #] Numerical factor :
1714*/
1715 else {
1716oneterm:;
1717 if ( factor == 0 || *factor > 2 ) {
1718 if ( pow > 0 ) {
1719 *r1++ = -SNUMBER; *r1++ = -1;
1720 t = r5;
1721 while ( t < r3 ) {
1722 t += *t; t[-1] = -t[-1];
1723 }
1724 }
1725 t = rr; *r1++ = *t++; *r1++ = 1; t++;
1726 COPYARG(r1,t);
1727 while ( t < m ) *r1++ = *t++;
1728 }
1729 }
1730 TermFree(EAscrat,"execarg");
1731 }
1732 } /* AC.OldFactArgFlag */
1733 }
1734/* r1 is fout in ons voorbeeld. */
1735 r2[1] = r1 - r2;
1736 action = 1;
1737 v[2] = DIRTYFLAG;
1738 }
1739 else {
1740 r = t + t[1];
1741 while ( t < r ) *r1++ = *t++;
1742 }
1743 }
1744 r = term + *term;
1745 while ( t < r ) *r1++ = *t++;
1746 m = AT.WorkPointer;
1747 i = m[0] = r1 - m;
1748 t = term;
1749 while ( --i >= 0 ) *t++ = *m++;
1750 if ( AT.WorkPointer < t ) AT.WorkPointer = t;
1751 }
1752/*
1753 #] FACTARG :
1754*/
1755 AR.Cnumlhs = oldnumlhs;
1756 if ( action && Normalize(BHEAD term) ) goto execargerr;
1757 AT.WorkPointer = oldwork;
1758 if ( AT.WorkPointer < term + *term ) AT.WorkPointer = term + *term;
1759 AT.pWorkPointer = oldppointer;
1760 return(action);
1761execargerr:
1762 AT.WorkPointer = oldwork;
1763 AT.pWorkPointer = oldppointer;
1764 MLOCK(ErrorMessageLock);
1765 MesCall("execarg");
1766 MUNLOCK(ErrorMessageLock);
1767 return(-1);
1768}
1769
1770/*
1771 #] execarg :
1772 #[ execterm :
1773*/
1774
1775int execterm(PHEAD WORD *term, WORD level)
1776{
1777 GETBIDENTITY
1778 CBUF *C = cbuf+AM.rbufnum;
1779 WORD oldnumlhs = AR.Cnumlhs;
1780 WORD maxisat = C->lhs[level][2];
1781 WORD *buffer1 = 0;
1782 WORD *oldworkpointer = AT.WorkPointer;
1783 WORD *t1, i;
1784 WORD olddeferflag = AR.DeferFlag, tryterm = 0;
1785 AR.DeferFlag = 0;
1786 do {
1787 AR.Cnumlhs = C->lhs[level][3];
1788 NewSort(BHEAD0);
1789/*
1790 Normally for function arguments we do not use PolyFun/PolyRatFun.
1791 Hence NewSort sets the corresponding variables to zero.
1792 Here we overwrite that.
1793*/
1794 AN.FunSorts[AR.sLevel]->PolyFlag = ( AR.PolyFun != 0 ) ? AR.PolyFunType: 0;
1795 if ( AR.PolyFun == 0 ) { AN.FunSorts[AR.sLevel]->PolyFlag = 0; }
1796 else if ( AR.PolyFunType == 1 ) { AN.FunSorts[AR.sLevel]->PolyFlag = 1; }
1797 else if ( AR.PolyFunType == 2 ) {
1798 if ( AR.PolyFunExp == 2 ) AN.FunSorts[AR.sLevel]->PolyFlag = 1;
1799 else AN.FunSorts[AR.sLevel]->PolyFlag = 2;
1800 }
1801 if ( buffer1 ) {
1802 term = buffer1;
1803 while ( *term ) {
1804 t1 = oldworkpointer;
1805 i = *term; while ( --i >= 0 ) *t1++ = *term++;
1806 AT.WorkPointer = t1;
1807 if ( Generator(BHEAD oldworkpointer,level) ) goto exectermerr;
1808 }
1809 }
1810 else {
1811 if ( Generator(BHEAD term,level) ) goto exectermerr;
1812 }
1813 if ( buffer1 ) {
1814 if ( tryterm ) { TermFree(buffer1,"buffer in sort statement"); tryterm = 0; }
1815 else { M_free((void *)buffer1,"buffer in sort statement"); }
1816 buffer1 = 0;
1817 }
1818 AN.tryterm = 1;
1819 if ( EndSort(BHEAD (WORD *)((void *)(&buffer1)),2) < 0 ) goto exectermerr;
1820 tryterm = AN.tryterm; AN.tryterm = 0;
1821 level = AR.Cnumlhs;
1822 } while ( AR.Cnumlhs < maxisat );
1823 AR.Cnumlhs = oldnumlhs;
1824 AR.DeferFlag = olddeferflag;
1825 term = buffer1;
1826 while ( *term ) {
1827 t1 = oldworkpointer;
1828 i = *term; while ( --i >= 0 ) *t1++ = *term++;
1829 AT.WorkPointer = t1;
1830 if ( Generator(BHEAD oldworkpointer,level) ) goto exectermerr;
1831 }
1832 if ( tryterm ) { TermFree(buffer1,"buffer in term statement"); tryterm = 0; }
1833 else { M_free(buffer1,"buffer in term statement"); }
1834 buffer1 = 0;
1835 AT.WorkPointer = oldworkpointer;
1836 return(0);
1837exectermerr:
1838 AT.WorkPointer = oldworkpointer;
1839 AR.DeferFlag = olddeferflag;
1840 MLOCK(ErrorMessageLock);
1841 MesCall("execterm");
1842 MUNLOCK(ErrorMessageLock);
1843 return(-1);
1844}
1845
1846/*
1847 #] execterm :
1848 #[ ArgumentImplode :
1849*/
1850
1851int ArgumentImplode(PHEAD WORD *term, WORD *thelist)
1852{
1853 GETBIDENTITY
1854 WORD *liststart, *liststop, *inlist;
1855 WORD *w, *t, *tend, *tstop, *tt, *ttstop, *ttt, ncount, i;
1856 int action = 0;
1857 liststop = thelist + thelist[1];
1858 liststart = thelist + 2;
1859 t = term;
1860 tend = t + *t;
1861 tstop = tend - ABS(tend[-1]);
1862 t++;
1863 while ( t < tstop ) {
1864 if ( *t >= FUNCTION ) {
1865 inlist = liststart;
1866 while ( inlist < liststop && *inlist != *t ) inlist += inlist[1];
1867 if ( inlist < liststop ) {
1868 tt = t; ttstop = t + t[1]; w = AT.WorkPointer;
1869 for ( i = 0; i < FUNHEAD; i++ ) *w++ = *tt++;
1870 while ( tt < ttstop ) {
1871 ncount = 0;
1872 if ( *tt == -SNUMBER && tt[1] == 0 ) {
1873 ncount = 1; ttt = tt; tt += 2;
1874 while ( tt < ttstop && *tt == -SNUMBER && tt[1] == 0 ) {
1875 ncount++; tt += 2;
1876 }
1877 }
1878 if ( ncount > 0 ) {
1879 if ( tt < ttstop && *tt == -SNUMBER && ( tt[1] == 1 || tt[1] == -1 ) ) {
1880 *w++ = -SNUMBER;
1881 *w++ = (ncount+1) * tt[1];
1882 tt += 2;
1883 action = 1;
1884 }
1885 else if ( ( tt[0] == tt[ARGHEAD] + ARGHEAD )
1886 && ( ABS(tt[tt[0]-1]) == 3 )
1887 && ( tt[tt[0]-2] == 1 )
1888 && ( tt[tt[0]-3] == 1 ) ) { /* Single term with coef +/- 1 */
1889 i = *tt; NCOPY(w,tt,i)
1890 w[-3] = ncount+1;
1891 action = 1;
1892 }
1893 else if ( *tt == -SYMBOL ) {
1894 *w++ = ARGHEAD+8;
1895 *w++ = 0;
1896 FILLARG(w)
1897 *w++ = 8;
1898 *w++ = SYMBOL;
1899 *w++ = tt[1];
1900 *w++ = 1;
1901 *w++ = ncount+1; *w++ = 1; *w++ = 3;
1902 tt += 2;
1903 action = 1;
1904 }
1905 else if ( *tt <= -FUNCTION ) {
1906 *w++ = ARGHEAD+FUNHEAD+4;
1907 *w++ = 0;
1908 FILLARG(w)
1909 *w++ = -*tt++;
1910 *w++ = FUNHEAD+4;
1911 FILLFUN(w)
1912 *w++ = ncount+1; *w++ = 1; *w++ = 3;
1913 action = 1;
1914 }
1915 else {
1916 while ( ttt < tt ) *w++ = *ttt++;
1917 if ( tt < ttstop && *tt == -SNUMBER ) {
1918 *w++ = *tt++; *w++ = *tt++;
1919 }
1920 }
1921 }
1922 else if ( *tt <= -FUNCTION ) {
1923 *w++ = *tt++;
1924 }
1925 else if ( *tt < 0 ) {
1926 *w++ = *tt++;
1927 *w++ = *tt++;
1928 }
1929 else {
1930 i = *tt; NCOPY(w,tt,i)
1931 }
1932 }
1933 AT.WorkPointer[1] = w - AT.WorkPointer;
1934 while ( tt < tend ) *w++ = *tt++;
1935 ttt = AT.WorkPointer; tt = t;
1936 while ( ttt < w ) *tt++ = *ttt++;
1937 term[0] = tt - term;
1938 AT.WorkPointer = tt;
1939 tend = tt; tstop = tt - ABS(tt[-1]);
1940 }
1941 }
1942 t += t[1];
1943 }
1944 if ( action ) {
1945 if ( Normalize(BHEAD term) ) return(-1);
1946 }
1947 return(0);
1948}
1949
1950/*
1951 #] ArgumentImplode :
1952 #[ ArgumentExplode :
1953*/
1954
1955int ArgumentExplode(PHEAD WORD *term, WORD *thelist)
1956{
1957 GETBIDENTITY
1958 WORD *liststart, *liststop, *inlist, *old;
1959 WORD *w, *t, *tend, *tstop, *tt, *ttstop, *ttt, ncount, i;
1960 int action = 0;
1961 LONG x;
1962 liststop = thelist + thelist[1];
1963 liststart = thelist + 2;
1964 t = term;
1965 tend = t + *t;
1966 tstop = tend - ABS(tend[-1]);
1967 t++;
1968 while ( t < tstop ) {
1969 if ( *t >= FUNCTION ) {
1970 inlist = liststart;
1971 while ( inlist < liststop && *inlist != *t ) inlist += inlist[1];
1972 if ( inlist < liststop ) {
1973 tt = t; ttstop = t + t[1]; w = AT.WorkPointer;
1974 for ( i = 0; i < FUNHEAD; i++ ) *w++ = *tt++;
1975 while ( tt < ttstop ) {
1976 if ( *tt == -SNUMBER && tt[1] != 0 ) {
1977 if ( tt[1] < AM.MaxTer/((WORD)sizeof(WORD)*4)
1978 && tt[1] > -(AM.MaxTer/((WORD)sizeof(WORD)*4))
1979 && ( tt[1] > 1 || tt[1] < -1 ) ) {
1980 ncount = ABS(tt[1]);
1981 while ( ncount > 1 ) {
1982 *w++ = -SNUMBER; *w++ = 0; ncount--;
1983 }
1984 *w++ = -SNUMBER;
1985 if ( tt[1] < 0 ) *w++ = -1;
1986 else *w++ = 1;
1987 tt += 2;
1988 action = 1;
1989 }
1990 else {
1991 *w++ = *tt++; *w++ = *tt++;
1992 }
1993 }
1994 else if ( *tt <= -FUNCTION ) {
1995 *w++ = *tt++;
1996 }
1997 else if ( *tt < 0 ) {
1998 *w++ = *tt++;
1999 *w++ = *tt++;
2000 }
2001 else if ( tt[0] == tt[ARGHEAD]+ARGHEAD ) {
2002 ttt = tt + tt[0] - 1;
2003 i = (ABS(ttt[0])-1)/2;
2004 if ( i > 1 ) {
2005TooMany: old = AN.currentTerm;
2006 AN.currentTerm = term;
2007 MesPrint("Too many arguments in output of ArgExplode");
2008 MesPrint("Term = %t");
2009 AN.currentTerm = old;
2010 return(-1);
2011 }
2012 if ( ttt[-1] != 1 ) goto NoExplode;
2013 x = ttt[-2];
2014 if ( 2*x > (AT.WorkTop-w)-*term ) goto TooMany;
2015 ncount = x - 1;
2016 while ( ncount > 0 ) {
2017 *w++ = -SNUMBER; *w++ = 0; ncount--;
2018 }
2019 ttt[-2] = 1;
2020 i = *tt; NCOPY(w,tt,i)
2021 action = 1;
2022 }
2023 else {
2024NoExplode:
2025 i = *tt; NCOPY(w,tt,i)
2026 }
2027 }
2028 AT.WorkPointer[1] = w - AT.WorkPointer;
2029 while ( tt < tend ) *w++ = *tt++;
2030 ttt = AT.WorkPointer; tt = t;
2031 while ( ttt < w ) *tt++ = *ttt++;
2032 term[0] = tt - term;
2033 AT.WorkPointer = tt;
2034 tend = tt; tstop = tt - ABS(tt[-1]);
2035 }
2036 }
2037 t += t[1];
2038 }
2039 if ( action ) {
2040 if ( Normalize(BHEAD term) ) return(-1);
2041 }
2042 return(0);
2043}
2044
2045/*
2046 #] ArgumentExplode :
2047 #[ ArgFactorize :
2048*/
2064#define NEWORDER
2065
2066int ArgFactorize(PHEAD WORD *argin, WORD *argout)
2067{
2068/*
2069 #[ step 0 : Declarations and initializations
2070*/
2071 WORD *argfree, *argextra, *argcopy, *t, *tstop, *a, *a1, *a2;
2072#ifdef NEWORDER
2073 WORD *tt;
2074#endif
2075 WORD startebuf = cbuf[AT.ebufnum].numrhs,oldword;
2076 WORD oldsorttype = AR.SortType, numargs;
2077 int error = 0, action = 0, i, ii, number, sign = 1;
2078
2079 *argout = 0;
2080/*
2081 #] step 0 :
2082 #[ step 1 : Take care of ordering
2083*/
2084 AR.SortType = SORTHIGHFIRST;
2085 if ( oldsorttype != AR.SortType ) {
2086 NewSort(BHEAD0);
2087 oldword = argin[*argin]; argin[*argin] = 0;
2088 t = argin+ARGHEAD;
2089 while ( *t ) {
2090 tstop = t + *t;
2091 if ( AN.ncmod != 0 ) {
2092 if ( AN.ncmod != 1 || ( (WORD)AN.cmod[0] < 0 ) ) {
2093 MLOCK(ErrorMessageLock);
2094 MesPrint("Factorization modulus a number, greater than a WORD not implemented.");
2095 MUNLOCK(ErrorMessageLock);
2096 Terminate(-1);
2097 }
2098 if ( Modulus(t) ) {
2099 MLOCK(ErrorMessageLock);
2100 MesCall("ArgFactorize");
2101 MUNLOCK(ErrorMessageLock);
2102 Terminate(-1);
2103 }
2104 if ( !*t) { t = tstop; continue; }
2105 }
2106 StoreTerm(BHEAD t);
2107 t = tstop;
2108 }
2109 /* par = 1, in case the arg has more than SubTermsInSmall terms */
2110 EndSort(BHEAD argin+ARGHEAD,1);
2111 argin[*argin] = oldword;
2112 }
2113/*
2114 #] step 1 :
2115 #[ step 2 : take out the 'content'.
2116*/
2117 argfree = TakeArgContent(BHEAD argin,argout);
2118 {
2119 a1 = argout;
2120 while ( *a1 ) {
2121 if ( a1[0] == -SNUMBER && ( a1[1] == 1 || a1[1] == -1 ) ) {
2122 if ( a1[1] == -1 ) { sign = -sign; a1[1] = 1; }
2123 if ( a1[2] ) {
2124 a = t = a1+2; while ( *t ) NEXTARG(t);
2125 i = t - a1-2;
2126 t = a1; NCOPY(t,a,i);
2127 *t = 0;
2128 continue;
2129 }
2130 else {
2131 a1[0] = 0;
2132 }
2133 break;
2134 }
2135 else if ( a1[0] == FUNHEAD+ARGHEAD+4 && a1[ARGHEAD] == FUNHEAD+4
2136 && a1[*a1-1] == 3 && a1[*a1-2] == 1 && a1[*a1-3] == 1
2137 && a1[ARGHEAD+1] >= FUNCTION ) {
2138 a = t = a1+*a1; while ( *t ) NEXTARG(t);
2139 i = t - a;
2140 *a1 = -a1[ARGHEAD+1]; t = a1+1; NCOPY(t,a,i);
2141 *t = 0;
2142 }
2143 NEXTARG(a1);
2144 }
2145 }
2146 if ( argfree == 0 ) {
2147 argfree = argin;
2148 }
2149 else if ( argfree[0] == ( argfree[ARGHEAD]+ARGHEAD ) ) {
2150 Normalize(BHEAD argfree+ARGHEAD);
2151 argfree[0] = argfree[ARGHEAD]+ARGHEAD;
2152 argfree[1] = 0;
2153 if ( ( argfree[0] == ARGHEAD+4 ) && ( argfree[ARGHEAD+3] == 3 )
2154 && ( argfree[ARGHEAD+1] == 1 ) && ( argfree[ARGHEAD+2] == 1 ) ) {
2155 goto return0;
2156 }
2157 }
2158 else {
2159/*
2160 The way we took out objects is rather brutish. We have to
2161 normalize
2162*/
2163 NewSort(BHEAD0);
2164 t = argfree+ARGHEAD;
2165 while ( *t ) {
2166 tstop = t + *t;
2167 Normalize(BHEAD t);
2168 StoreTerm(BHEAD t);
2169 t = tstop;
2170 }
2171 /* par = 1, in case the arg has more than SubTermsInSmall terms */
2172 EndSort(BHEAD argfree+ARGHEAD,1);
2173 t = argfree+ARGHEAD;
2174 while ( *t ) t += *t;
2175 *argfree = t - argfree;
2176 }
2177/*
2178 #] step 2 :
2179 #[ step 3 : look whether we have done this one already.
2180*/
2181 if ( ( number = FindArg(BHEAD argfree) ) != 0 ) {
2182 if ( number > 0 ) t = cbuf[AT.fbufnum].rhs[number-1];
2183 else t = cbuf[AC.ffbufnum].rhs[-number-1];
2184/*
2185 Now position on the result. Remember we have in the cache:
2186 inputarg,0,outputargs,0
2187 t is currently at inputarg. *inputarg is always positive.
2188 in principle this holds also for the arguments in the output
2189 but we take no risks here (in case of future developments).
2190*/
2191 t += *t; t++;
2192 tstop = t;
2193 ii = 0;
2194 while ( *tstop ) {
2195 if ( *tstop == -SNUMBER && tstop[1] == -1 ) {
2196 sign = -sign; ii += 2;
2197 }
2198 NEXTARG(tstop);
2199 }
2200 a = argout; while ( *a ) NEXTARG(a);
2201#ifndef NEWORDER
2202 if ( sign == -1 ) { *a++ = -SNUMBER; *a++ = -1; *a = 0; sign = 1; }
2203#endif
2204 i = tstop - t - ii;
2205 ii = a - argout;
2206 a2 = a; a1 = a + i;
2207 *a1 = 0;
2208 while ( ii > 0 ) { *--a1 = *--a2; ii--; }
2209 a = argout;
2210 while ( *t ) {
2211 if ( *t == -SNUMBER && t[1] == -1 ) { t += 2; }
2212 else { COPY1ARG(a,t) }
2213 }
2214 goto return0;
2215 }
2216/*
2217 #] step 3 :
2218 #[ step 4 : invoke ConvertToPoly
2219
2220 We make a copy first in case there are no factors
2221*/
2222 argcopy = TermMalloc("argcopy");
2223 for ( i = 0; i <= *argfree; i++ ) argcopy[i] = argfree[i];
2224
2225 tstop = argfree + *argfree;
2226 {
2227 WORD sumcommu = 0;
2228 t = argfree + ARGHEAD;
2229 while ( t < tstop ) {
2230 sumcommu += DoesCommu(t);
2231 t += *t;
2232 }
2233 if ( sumcommu > 1 ) {
2234 MesPrint("ERROR: Cannot factorize an argument with more than one noncommuting object");
2235 Terminate(-1);
2236 }
2237 }
2238 t = argfree + ARGHEAD;
2239
2240 while ( t < tstop ) {
2241 if ( ( t[1] != SYMBOL ) && ( *t != (ABS(t[*t-1])+1) ) ) {
2242 action = 1; break;
2243 }
2244 t += *t;
2245 }
2246 if ( action ) {
2247 t = argfree + ARGHEAD;
2248 argextra = AT.WorkPointer;
2249 NewSort(BHEAD0);
2250 while ( t < tstop ) {
2251 if ( LocalConvertToPoly(BHEAD t,argextra,startebuf,0) < 0 ) {
2252 error = -1;
2253getout:
2254 AR.SortType = oldsorttype;
2255 TermFree(argcopy,"argcopy");
2256 if ( argfree != argin ) TermFree(argfree,"argfree");
2257 MesCall("ArgFactorize");
2258 Terminate(-1);
2259 return(-1);
2260 }
2261 StoreTerm(BHEAD argextra);
2262 t += *t; argextra += *argextra;
2263 }
2264 /* par = 1, in case the arg has more than SubTermsInSmall terms */
2265 if ( EndSort(BHEAD argfree+ARGHEAD,1) < 0 ) { error = -2; goto getout; }
2266 t = argfree + ARGHEAD;
2267 while ( *t > 0 ) t += *t;
2268 *argfree = t - argfree;
2269 }
2270/*
2271 #] step 4 :
2272 #[ step 5 : If not in the tables, we have to do this by hard work.
2273*/
2274
2275 a = argout;
2276 while ( *a ) NEXTARG(a);
2277 if ( poly_factorize_argument(BHEAD argfree,a) < 0 ) {
2278 MesCall("ArgFactorize");
2279 error = -1;
2280 }
2281/*
2282 #] step 5 :
2283 #[ step 6 : use now ConvertFromPoly
2284
2285 Be careful: there should be more than one argument now.
2286*/
2287 if ( error == 0 && action ) {
2288 a1 = a; NEXTARG(a1);
2289 if ( *a1 != 0 ) {
2290 CBUF *C = cbuf+AC.cbufnum;
2291 CBUF *CC = cbuf+AT.ebufnum;
2292 WORD *oldworkpointer = AT.WorkPointer;
2293 WORD *argcopy2 = TermMalloc("argcopy2"), *a1, *a2;
2294 a1 = a; a2 = argcopy2;
2295 while ( *a1 ) {
2296 if ( *a1 < 0 ) {
2297 if ( *a1 > -FUNCTION ) *a2++ = *a1++;
2298 *a2++ = *a1++; *a2 = 0;
2299 continue;
2300 }
2301 t = a1 + ARGHEAD;
2302 tstop = a1 + *a1;
2303 argextra = AT.WorkPointer;
2304 NewSort(BHEAD0);
2305 while ( t < tstop ) {
2306 if ( ConvertFromPoly(BHEAD t,argextra,numxsymbol,CC->numrhs-startebuf+numxsymbol
2307 ,startebuf-numxsymbol,1) <= 0 ) {
2308 TermFree(argcopy2,"argcopy2");
2309 LowerSortLevel();
2310 error = -3;
2311 goto getout;
2312 }
2313 t += *t;
2314 AT.WorkPointer = argextra + *argextra;
2315/*
2316 ConvertFromPoly leaves terms with subexpressions. Hence:
2317*/
2318 if ( Generator(BHEAD argextra,C->numlhs) ) {
2319 TermFree(argcopy2,"argcopy2");
2320 LowerSortLevel();
2321 error = -4;
2322 goto getout;
2323 }
2324 }
2325 AT.WorkPointer = oldworkpointer;
2326 /* par = 1, in case the factor has more than SubTermsInSmall terms */
2327 if ( EndSort(BHEAD a2+ARGHEAD,1) < 0 ) { error = -5; goto getout; }
2328 t = a2+ARGHEAD; while ( *t ) t += *t;
2329 *a2 = t - a2; a2[1] = 0; ZEROARG(a2);
2330 ToFast(a2,a2); NEXTARG(a2);
2331 a1 = tstop;
2332 }
2333 i = a2 - argcopy2;
2334 a2 = argcopy2; a1 = a;
2335 NCOPY(a1,a2,i);
2336 *a1 = 0;
2337 TermFree(argcopy2,"argcopy2");
2338/*
2339 Erase the entries we made temporarily in cbuf[AT.ebufnum]
2340*/
2341 CC->numrhs = startebuf;
2342 }
2343 else { /* no factorization. recover the argument from before step 3. */
2344 for ( i = 0; i <= *argcopy; i++ ) a[i] = argcopy[i];
2345 }
2346 }
2347/*
2348 #] step 6 :
2349 #[ step 7 : Add this one to the tables.
2350
2351 Possibly drop some elements in the tables
2352 when they become too full.
2353*/
2354 if ( error == 0 && AN.ncmod == 0 ) {
2355 if ( InsertArg(BHEAD argcopy,a,0) < 0 ) { error = -1; }
2356 }
2357/*
2358 #] step 7 :
2359 #[ step 8 : Clean up and return.
2360
2361 Change the order of the arguments in argout and a.
2362 Use argcopy as spare space.
2363*/
2364 ii = a - argout;
2365 for ( i = 0; i < ii; i++ ) argcopy[i] = argout[i];
2366 a1 = a;
2367 while ( *a1 ) {
2368 if ( *a1 == -SNUMBER && a1[1] < 0 ) {
2369 sign = -sign; a1[1] = -a1[1];
2370 if ( a1[1] == 1 ) {
2371 a2 = a1+2; while ( *a2 ) NEXTARG(a2);
2372 i = a2-a1-2; a2 = a1+2;
2373 NCOPY(a1,a2,i);
2374 *a1 = 0;
2375 }
2376 while ( *a1 ) NEXTARG(a1);
2377 break;
2378 }
2379 else {
2380 if ( *a1 > 0 && *a1 == a1[ARGHEAD]+ARGHEAD && a1[*a1-1] < 0 ) {
2381 a1[*a1-1] = -a1[*a1-1]; sign = -sign;
2382 }
2383 if ( *a1 == ARGHEAD+4 && a1[ARGHEAD+1] == 1 && a1[ARGHEAD+2] == 1 ) {
2384 a2 = a1+ARGHEAD+4; while ( *a2 ) NEXTARG(a2);
2385 i = a2-a1-ARGHEAD-4; a2 = a1+ARGHEAD+4;
2386 NCOPY(a1,a2,i);
2387 *a1 = 0;
2388 break;
2389 }
2390 while ( *a1 ) NEXTARG(a1);
2391 break;
2392 }
2393 NEXTARG(a1);
2394 }
2395 i = a1 - a;
2396 a2 = argout;
2397 NCOPY(a2,a,i);
2398 for ( i = 0; i < ii; i++ ) *a2++ = argcopy[i];
2399#ifndef NEWORDER
2400 if ( sign == -1 ) { *a2++ = -SNUMBER; *a2++ = -1; sign = 1; }
2401#endif
2402 *a2 = 0;
2403 TermFree(argcopy,"argcopy");
2404return0:
2405 if ( argfree != argin ) TermFree(argfree,"argfree");
2406 if ( oldsorttype != AR.SortType ) {
2407 AR.SortType = oldsorttype;
2408 a = argout;
2409 while ( *a ) {
2410 if ( *a > 0 ) {
2411 NewSort(BHEAD0);
2412 oldword = a[*a]; a[*a] = 0;
2413 t = a+ARGHEAD;
2414 while ( *t ) {
2415 tstop = t + *t;
2416 StoreTerm(BHEAD t);
2417 t = tstop;
2418 }
2419 /* par = 1, in case the factor has more than SubTermsInSmall terms */
2420 EndSort(BHEAD a+ARGHEAD,1);
2421 a[*a] = oldword;
2422 a += *a;
2423 }
2424 else { NEXTARG(a); }
2425 }
2426 }
2427#ifdef NEWORDER
2428 t = argout; numargs = 0;
2429 while ( *t ) {
2430 tt = t;
2431 NEXTARG(t);
2432 if ( *tt == ABS(t[-1])+1+ARGHEAD && sign == -1 ) { t[-1] = -t[-1]; sign = 1; }
2433 else if ( *tt == -SNUMBER && sign == -1 ) { tt[1] = -tt[1]; sign = 1; }
2434 numargs++;
2435 }
2436 if ( sign == -1 ) {
2437 *t++ = -SNUMBER; *t++ = -1; *t = 0; sign = 1; numargs++;
2438 }
2439#else
2440/*
2441 Now we have to sort the arguments
2442 First have the number of 'nontrivial/nonnumerical' arguments
2443 Then make a piece of code like in FullSymmetrize with that number
2444 of arguments to be symmetrized.
2445 Put a function in front
2446 Call the Symmetrize routine
2447*/
2448 t = argout; numargs = 0;
2449 while ( *t && *t != -SNUMBER && ( *t < 0 || ( ABS(t[*t-1]) != *t-1 ) ) ) {
2450 NEXTARG(t);
2451 numargs++;
2452 }
2453#endif
2454 if ( numargs > 1 ) {
2455 WORD *Lijst;
2456 WORD x[3];
2457 x[0] = argout[-FUNHEAD];
2458 x[1] = argout[-FUNHEAD+1];
2459 x[2] = argout[-FUNHEAD+2];
2460 while ( *t ) { NEXTARG(t); }
2461 argout[-FUNHEAD] = SQRTFUNCTION;
2462 argout[-FUNHEAD+1] = t-argout+FUNHEAD;
2463 argout[-FUNHEAD+2] = 0;
2464 AT.WorkPointer = t+1;
2465 Lijst = AT.WorkPointer;
2466 for ( i = 0; i < numargs; i++ ) Lijst[i] = i;
2467 AT.WorkPointer += numargs;
2468 error = Symmetrize(BHEAD argout-FUNHEAD,Lijst,numargs,1,SYMMETRIC);
2469 AT.WorkPointer = Lijst;
2470 argout[-FUNHEAD] = x[0];
2471 argout[-FUNHEAD+1] = x[1];
2472 argout[-FUNHEAD+2] = x[2];
2473#ifdef NEWORDER
2474/*
2475 Now we have to get a potential numerical argument to the first position
2476*/
2477 tstop = argout; while ( *tstop ) { NEXTARG(tstop); }
2478 t = argout; number = 0;
2479 while ( *t ) {
2480 tt = t; NEXTARG(t);
2481 if ( *tt == -SNUMBER ) {
2482 if ( number == 0 ) break;
2483 x[0] = tt[1];
2484 while ( tt > argout ) { *--t = *--tt; }
2485 argout[0] = -SNUMBER; argout[1] = x[0];
2486 break;
2487 }
2488 else if ( *tt == ABS(t[-1])+1+ARGHEAD ) {
2489 if ( number == 0 ) break;
2490 ii = t - tt;
2491 for ( i = 0; i < ii; i++ ) tstop[i] = tt[i];
2492 while ( tt > argout ) { *--t = *--tt; }
2493 for ( i = 0; i < ii; i++ ) argout[i] = tstop[i];
2494 *tstop = 0;
2495 break;
2496 }
2497 number++;
2498 }
2499#endif
2500 }
2501/*
2502 #] step 8 :
2503*/
2504 return(error);
2505}
2506
2507/*
2508 #] ArgFactorize :
2509 #[ FindArg :
2510*/
2519WORD FindArg(PHEAD WORD *a)
2520{
2521 int number;
2522 if ( AN.ncmod != 0 ) return(0); /* no room for mod stuff */
2523 number = FindTree(AT.fbufnum,a);
2524 if ( number >= 0 ) return(number+1);
2525 number = FindTree(AC.ffbufnum,a);
2526 if ( number >= 0 ) return(-number-1);
2527 return(0);
2528}
2529
2530/*
2531 #] FindArg :
2532 #[ InsertArg :
2533*/
2543int InsertArg(PHEAD WORD *argin, WORD *argout,int par)
2544{
2545 CBUF *C;
2546 WORD *a, i, bufnum;
2547 if ( par == 0 ) {
2548 bufnum = AT.fbufnum;
2549 C = cbuf+bufnum;
2550 if ( C->numrhs >= (C->maxrhs-2) ) CleanupArgCache(BHEAD AT.fbufnum);
2551 }
2552 else if ( par == 1 ) {
2553 bufnum = AC.ffbufnum;
2554 C = cbuf+bufnum;
2555 }
2556 else { return(-1); }
2557 AddRHS(bufnum,1);
2558 AddNtoC(bufnum,*argin,argin,1);
2559 AddToCB(C,0)
2560 a = argout; while ( *a ) NEXTARG(a);
2561 i = a - argout;
2562 AddNtoC(bufnum,i,argout,2);
2563 AddToCB(C,0)
2564 return(InsTree(bufnum,C->numrhs));
2565}
2566
2567/*
2568 #] InsertArg :
2569 #[ CleanupArgCache :
2570*/
2578int CleanupArgCache(PHEAD WORD bufnum)
2579{
2580 CBUF *C = cbuf + bufnum;
2581 COMPTREE *boomlijst = C->boomlijst;
2582 LONG *weights = (LONG *)Malloc1(2*(C->numrhs+1)*sizeof(LONG),"CleanupArgCache");
2583 LONG w, whalf, *extraweights;
2584 WORD *a, *to, *from;
2585 int i,j,k;
2586 for ( i = 1; i <= C->numrhs; i++ ) {
2587 weights[i] = ((LONG)i) * (boomlijst[i].usage);
2588 }
2589/*
2590 Now sort the weights and determine the halfway weight
2591*/
2592 extraweights = weights+C->numrhs+1;
2593 SortWeights(weights+1,extraweights,C->numrhs);
2594 whalf = weights[C->numrhs/2+1];
2595/*
2596 We should drop everybody with a weight < whalf.
2597*/
2598 to = C->Buffer;
2599 k = 1;
2600 for ( i = 1; i <= C->numrhs; i++ ) {
2601 from = C->rhs[i]; w = ((LONG)i) * (boomlijst[i].usage);
2602 if ( w >= whalf ) {
2603 if ( i < C->numrhs-1 ) {
2604 if ( to == from ) {
2605 to = C->rhs[i+1];
2606 }
2607 else {
2608 j = C->rhs[i+1] - from;
2609 C->rhs[k] = to;
2610 NCOPY(to,from,j)
2611 }
2612 }
2613 else if ( to == from ) {
2614 to += *to + 1; while ( *to ) NEXTARG(to); to++;
2615 }
2616 else {
2617 a = from; a += *a+1; while ( *a ) NEXTARG(a); a++;
2618 j = a - from;
2619 C->rhs[k] = to;
2620 NCOPY(to,from,j)
2621 }
2622 weights[k++] = boomlijst[i].usage;
2623 }
2624 }
2625 C->numrhs = --k;
2626 C->Pointer = to;
2627/*
2628 Next we need to rebuild the tree.
2629 Note that this can probably be done much faster by using the
2630 remains of the old tree !!!!!!!!!!!!!!!!
2631*/
2632 ClearTree(AT.fbufnum);
2633 for ( i = 1; i <= k; i++ ) {
2634 InsTree(AT.fbufnum,i);
2635 boomlijst[i].usage = weights[i];
2636 }
2637/*
2638 And cleanup
2639*/
2640 M_free(weights,"CleanupArgCache");
2641 return(0);
2642}
2643
2644/*
2645 #] CleanupArgCache :
2646 #[ ArgSymbolMerge :
2647*/
2648
2649int ArgSymbolMerge(WORD *t1, WORD *t2)
2650{
2651 WORD *t1e = t1+t1[1];
2652 WORD *t2e = t2+t2[1];
2653 WORD *t1a = t1+2;
2654 WORD *t2a = t2+2;
2655 WORD *t3;
2656 while ( t1a < t1e && t2a < t2e ) {
2657 if ( *t1a < *t2a ) {
2658 if ( t1a[1] >= 0 ) {
2659 t3 = t1a+2;
2660 while ( t3 < t1e ) { t3[-2] = *t3; t3[-1] = t3[1]; t3 += 2; }
2661 t1e -= 2;
2662 }
2663 else t1a += 2;
2664 }
2665 else if ( *t1a > *t2a ) {
2666 if ( t2a[1] >= 0 ) t2a += 2;
2667 else {
2668 t3 = t1e;
2669 while ( t3 > t1a ) { *t3 = t3[-2]; t3[1] = t3[-1]; t3 -= 2; }
2670 *t1a++ = *t2a++;
2671 *t1a++ = *t2a++;
2672 t1e += 2;
2673 }
2674 }
2675 else {
2676 if ( t2a[1] < t1a[1] ) t1a[1] = t2a[1];
2677 t1a += 2; t2a += 2;
2678 }
2679 }
2680 while ( t2a < t2e ) {
2681 if ( t2a[1] < 0 ) {
2682 *t1a++ = *t2a++;
2683 *t1a++ = *t2a++;
2684 }
2685 else t2a += 2;
2686 }
2687 while ( t1a < t1e ) {
2688 if ( t1a[1] >= 0 ) {
2689 t3 = t1a+2;
2690 while ( t3 < t1e ) { t3[-2] = *t3; t3[-1] = t3[1]; t3 += 2; }
2691 t1e -= 2;
2692 }
2693 else t1a += 2;
2694 }
2695 t1[1] = t1a - t1;
2696 return(0);
2697}
2698
2699/*
2700 #] ArgSymbolMerge :
2701 #[ ArgDotproductMerge :
2702*/
2703
2704int ArgDotproductMerge(WORD *t1, WORD *t2)
2705{
2706 WORD *t1e = t1+t1[1];
2707 WORD *t2e = t2+t2[1];
2708 WORD *t1a = t1+2;
2709 WORD *t2a = t2+2;
2710 WORD *t3;
2711 while ( t1a < t1e && t2a < t2e ) {
2712 if ( *t1a < *t2a || ( *t1a == *t2a && t1a[1] < t2a[1] ) ) {
2713 if ( t1a[2] >= 0 ) {
2714 t3 = t1a+3;
2715 while ( t3 < t1e ) { t3[-3] = *t3; t3[-2] = t3[1]; t3[-1] = t3[2]; t3 += 3; }
2716 t1e -= 3;
2717 }
2718 else t1a += 3;
2719 }
2720 else if ( *t1a > *t2a || ( *t1a == *t2a && t1a[1] > t2a[1] ) ) {
2721 if ( t2a[2] >= 0 ) t2a += 3;
2722 else {
2723 t3 = t1e;
2724 while ( t3 > t1a ) { *t3 = t3[-3]; t3[1] = t3[-2]; t3[2] = t3[-1]; t3 -= 3; }
2725 *t1a++ = *t2a++;
2726 *t1a++ = *t2a++;
2727 *t1a++ = *t2a++;
2728 t1e += 3;
2729 }
2730 }
2731 else {
2732 if ( t2a[2] < t1a[2] ) t1a[2] = t2a[2];
2733 t1a += 3; t2a += 3;
2734 }
2735 }
2736 while ( t2a < t2e ) {
2737 if ( t2a[2] < 0 ) {
2738 *t1a++ = *t2a++;
2739 *t1a++ = *t2a++;
2740 *t1a++ = *t2a++;
2741 }
2742 else t2a += 3;
2743 }
2744 while ( t1a < t1e ) {
2745 if ( t1a[2] >= 0 ) {
2746 t3 = t1a+3;
2747 while ( t3 < t1e ) { t3[-3] = *t3; t3[-2] = t3[1]; t3[-1] = t3[2]; t3 += 3; }
2748 t1e -= 3;
2749 }
2750 else t1a += 2;
2751 }
2752 t1[1] = t1a - t1;
2753 return(0);
2754}
2755
2756/*
2757 #] ArgDotproductMerge :
2758 #[ TakeArgContent :
2759*/
2772WORD *TakeArgContent(PHEAD WORD *argin, WORD *argout)
2773{
2774 GETBIDENTITY
2775 WORD *t, *rnext, *r1, *r2, *r3, *r5, *r6, *r7, *r8, *r9;
2776 WORD pow, *mm, *mnext, *mstop, *argin2 = argin, *argin3 = argin, *argfree;
2777 WORD ncom;
2778 int j, i, act;
2779 r5 = t = argin + ARGHEAD;
2780 r3 = argin + *argin;
2781 rnext = t + *t;
2782 GETSTOP(t,r6);
2783 r1 = argout;
2784 t++;
2785/*
2786 First pass: arrange everything but the symbols and dotproducts.
2787 They need separate treatment because we have to take out negative
2788 powers.
2789*/
2790 while ( t < r6 ) {
2791/*
2792 #[ DELTA/VECTOR :
2793*/
2794 if ( *t == DELTA || *t == VECTOR ) {
2795 r7 = t; r8 = t + t[1]; t += 2;
2796 while ( t < r8 ) {
2797 mm = rnext;
2798 pow = 1;
2799 while ( mm < r3 ) {
2800 mnext = mm + *mm;
2801 GETSTOP(mm,mstop); mm++;
2802 while ( mm < mstop ) {
2803 if ( *mm != *r7 ) mm += mm[1];
2804 else break;
2805 }
2806 if ( *mm == *r7 ) {
2807 mstop = mm + mm[1]; mm += 2;
2808 while ( ( *mm != *t || mm[1] != t[1] )
2809 && mm < mstop ) mm += 2;
2810 if ( mm >= mstop ) pow = 0;
2811 }
2812 else pow = 0;
2813 if ( pow == 0 ) break;
2814 mm = mnext;
2815 }
2816 if ( pow == 0 ) { t += 2; continue; }
2817/*
2818 We have a factor
2819*/
2820 *r1++ = 8 + ARGHEAD;
2821 for ( j = 1; j < ARGHEAD; j++ ) *r1++ = 0;
2822 *r1++ = 8; *r1++ = *r7;
2823 *r1++ = 4; *r1++ = *t; *r1++ = t[1];
2824 *r1++ = 1; *r1++ = 1; *r1++ = 3;
2825 argout = r1;
2826/*
2827 Now we have to remove the delta's/vectors
2828*/
2829 mm = rnext;
2830 while ( mm < r3 ) {
2831 mnext = mm + *mm;
2832 GETSTOP(mm,mstop); mm++;
2833 while ( mm < mstop ) {
2834 if ( *mm != *r7 ) mm += mm[1];
2835 else break;
2836 }
2837 mstop = mm + mm[1]; mm += 2;
2838 while ( mm < mstop && (
2839 *mm != *t || mm[1] != t[1] ) ) mm += 2;
2840 *mm = mm[1] = NOINDEX;
2841 mm = mnext;
2842 }
2843 *t = t[1] = NOINDEX;
2844 t += 2;
2845 }
2846 }
2847/*
2848 #] DELTA/VECTOR :
2849 #[ INDEX :
2850*/
2851 else if ( *t == INDEX ) {
2852 r7 = t; r8 = t + t[1]; t += 2;
2853 while ( t < r8 ) {
2854 mm = rnext;
2855 pow = 1;
2856 while ( mm < r3 ) {
2857 mnext = mm + *mm;
2858 GETSTOP(mm,mstop); mm++;
2859 while ( mm < mstop ) {
2860 if ( *mm != *r7 ) mm += mm[1];
2861 else break;
2862 }
2863 if ( *mm == *r7 ) {
2864 mstop = mm + mm[1]; mm += 2;
2865 while ( *mm != *t
2866 && mm < mstop ) mm++;
2867 if ( mm >= mstop ) pow = 0;
2868 }
2869 else pow = 0;
2870 if ( pow == 0 ) break;
2871 mm = mnext;
2872 }
2873 if ( pow == 0 ) { t++; continue; }
2874/*
2875 We have a factor
2876*/
2877 if ( *t < 0 ) { *r1++ = -VECTOR; }
2878 else { *r1++ = -INDEX; }
2879 *r1++ = *t;
2880 argout = r1;
2881/*
2882 Now we have to remove the index
2883*/
2884 *t = NOINDEX;
2885 mm = rnext;
2886 while ( mm < r3 ) {
2887 mnext = mm + *mm;
2888 GETSTOP(mm,mstop); mm++;
2889 while ( mm < mstop ) {
2890 if ( *mm != *r7 ) mm += mm[1];
2891 else break;
2892 }
2893 mstop = mm + mm[1]; mm += 2;
2894 while ( mm < mstop &&
2895 *mm != *t ) mm += 1;
2896 *mm = NOINDEX;
2897 mm = mnext;
2898 }
2899 t += 1;
2900 }
2901 }
2902/*
2903 #] INDEX :
2904 #[ FUNCTION :
2905*/
2906 else if ( *t >= FUNCTION ) {
2907/*
2908 In the next code we should actually look inside
2909 the DENOMINATOR or EXPONENT for noncommuting objects
2910*/
2911 if ( *t >= FUNCTION &&
2912 functions[*t-FUNCTION].commute == 0 ) ncom = 0;
2913 else ncom = 1;
2914 if ( ncom ) {
2915 mm = r5 + 1;
2916 while ( mm < t && ( *mm == DUMMYFUN
2917 || *mm == DUMMYTEN ) ) mm += mm[1];
2918 if ( mm < t ) { t += t[1]; continue; }
2919 }
2920 mm = rnext; pow = 1;
2921 while ( mm < r3 ) {
2922 mnext = mm + *mm;
2923 GETSTOP(mm,mstop); mm++;
2924 while ( mm < mstop ) {
2925 if ( *mm == *t && mm[1] == t[1] ) {
2926 for ( i = 2; i < t[1]; i++ ) {
2927 if ( mm[i] != t[i] ) break;
2928 }
2929 if ( i >= t[1] )
2930 { mm += mm[1]; goto nextmterm; }
2931 }
2932 if ( ncom && *mm != DUMMYFUN && *mm != DUMMYTEN )
2933 { pow = 0; break; }
2934 mm += mm[1];
2935 }
2936 if ( mm >= mstop ) pow = 0;
2937 if ( pow == 0 ) break;
2938nextmterm: mm = mnext;
2939 }
2940 if ( pow == 0 ) { t += t[1]; continue; }
2941/*
2942 Copy the function
2943*/
2944 *r1++ = t[1] + 4 + ARGHEAD;
2945 for ( i = 1; i < ARGHEAD; i++ ) *r1++ = 0;
2946 *r1++ = t[1] + 4;
2947 for ( i = 0; i < t[1]; i++ ) *r1++ = t[i];
2948 *r1++ = 1; *r1++ = 1; *r1++ = 3;
2949 argout = r1;
2950/*
2951 Now we have to take out the functions
2952*/
2953 mm = rnext;
2954 while ( mm < r3 ) {
2955 mnext = mm + *mm;
2956 GETSTOP(mm,mstop); mm++;
2957 while ( mm < mstop ) {
2958 if ( *mm == *t && mm[1] == t[1] ) {
2959 for ( i = 2; i < t[1]; i++ ) {
2960 if ( mm[i] != t[i] ) break;
2961 }
2962 if ( i >= t[1] ) {
2963 if ( functions[*t-FUNCTION].spec > 0 )
2964 *mm = DUMMYTEN;
2965 else
2966 *mm = DUMMYFUN;
2967 mm += mm[1];
2968 goto nextterm;
2969 }
2970 }
2971 mm += mm[1];
2972 }
2973nextterm: mm = mnext;
2974 }
2975 if ( functions[*t-FUNCTION].spec > 0 )
2976 *t = DUMMYTEN;
2977 else
2978 *t = DUMMYFUN;
2979 t += t[1];
2980 }
2981/*
2982 #] FUNCTION :
2983*/
2984 else {
2985 t += t[1];
2986 }
2987 }
2988/*
2989 #[ SYMBOL :
2990
2991 Now collect all symbols. We can use the space after r1 as storage
2992*/
2993 t = argin+ARGHEAD;
2994 rnext = t + *t;
2995 r2 = r1;
2996 while ( t < r3 ) {
2997 GETSTOP(t,r6);
2998 t++;
2999 act = 0;
3000 while ( t < r6 ) {
3001 if ( *t == SYMBOL ) {
3002 act = 1;
3003 i = t[1];
3004 NCOPY(r2,t,i)
3005 }
3006 else { t += t[1]; }
3007 }
3008 if ( act == 0 ) {
3009 *r2++ = SYMBOL; *r2++ = 2;
3010 }
3011 t = rnext; rnext = rnext + *rnext;
3012 }
3013 *r2 = 0;
3014 argin2 = argin;
3015/*
3016 Now we have a list of all symbols as a sequence of SYMBOL subterms.
3017 Any symbol that is absent in a subterm has power zero.
3018 We now need a list of all minimum powers.
3019 This can be done by subsequent merges.
3020*/
3021 r7 = r1; /* The first object into which we merge. */
3022 r8 = r7 + r7[1]; /* The object that gets merged into r7. */
3023 while ( *r8 ) {
3024 r2 = r8 + r8[1]; /* Next object */
3025 ArgSymbolMerge(r7,r8);
3026 r8 = r2;
3027 }
3028/*
3029 Now we have to divide by the object in r7 and take it apart as factors.
3030 The division can be simple if there are no negative powers.
3031*/
3032 if ( r7[1] > 2 ) {
3033 r8 = r7+2;
3034 r2 = r7 + r7[1];
3035 act = 0;
3036 pow = 0;
3037 while ( r8 < r2 ) {
3038 if ( r8[1] < 0 ) { act = 1; pow += -r8[1]*(ARGHEAD+8); }
3039 else { pow += 2*r8[1]; }
3040 r8 += 2;
3041 }
3042/*
3043 The amount of space we need to move r7 is given in pow
3044*/
3045 if ( act == 0 ) { /* this can be done 'in situ' */
3046 t = argin + ARGHEAD;
3047 while ( t < r3 ) {
3048 rnext = t + *t;
3049 GETSTOP(t,r6);
3050 t++;
3051 while ( t < r6 ) {
3052 if ( *t != SYMBOL ) { t += t[1]; continue; }
3053 r8 = r7+2; r9 = t + t[1]; t += 2;
3054 while ( ( t < r9 ) && ( r8 < r2 ) ) {
3055 if ( *t == *r8 ) {
3056 t[1] -= r8[1]; t += 2; r8 += 2;
3057 }
3058 else { /* *t must be < than *r8 !!! */
3059 t += 2;
3060 }
3061 }
3062 t = r9;
3063 }
3064 t = rnext;
3065 }
3066/*
3067 And now the factors that go to argout.
3068 First we have to move r7 out of the way.
3069*/
3070 r8 = r7+pow; i = r7[1];
3071 while ( --i >= 0 ) r8[i] = r7[i];
3072 r2 += pow;
3073 r8 += 2;
3074 while ( r8 < r2 ) {
3075 for ( i = 0; i < r8[1]; i++ ) { *r1++ = -SYMBOL; *r1++ = *r8; }
3076 r8 += 2;
3077 }
3078 }
3079 else { /* this needs a new location */
3080 argin2 = TermMalloc("TakeArgContent2");
3081/*
3082 We have to multiply the inverse of r7 into argin
3083 The answer should go to argin2.
3084*/
3085 r5 = argin2; *r5++ = 0; *r5++ = 0; FILLARG(r5);
3086 t = argin+ARGHEAD;
3087 while ( t < r3 ) {
3088 rnext = t + *t;
3089 GETSTOP(t,r6);
3090 r9 = r5;
3091 *r5++ = *t++ + r7[1];
3092 while ( t < r6 ) *r5++ = *t++;
3093 i = r7[1] - 2; r8 = r7+2;
3094 *r5++ = r7[0]; *r5++ = r7[1];
3095 while ( i > 0 ) { *r5++ = *r8++; *r5++ = -*r8++; i -= 2; }
3096 while ( t < rnext ) *r5++ = *t++;
3097 Normalize(BHEAD r9);
3098 r5 = r9 + *r9;
3099 }
3100 *r5 = 0;
3101 *argin2 = r5-argin2;
3102/*
3103 We may have to sort the terms in argin2.
3104*/
3105 NewSort(BHEAD0);
3106 t = argin2+ARGHEAD;
3107 while ( *t ) {
3108 StoreTerm(BHEAD t);
3109 t += *t;
3110 }
3111 t = argin2+ARGHEAD;
3112 if ( EndSort(BHEAD t,0) < 0 ) goto Irreg;
3113 while ( *t ) t += *t;
3114 *argin2 = t - argin2;
3115 r3 = t;
3116/*
3117 And now the factors that go to argout.
3118 First we have to move r7 out of the way.
3119*/
3120 r8 = r7+pow; i = r7[1];
3121 while ( --i >= 0 ) r8[i] = r7[i];
3122 r2 += pow;
3123 r8 += 2;
3124 while ( r8 < r2 ) {
3125 if ( r8[1] >= 0 ) {
3126 for ( i = 0; i < r8[1]; i++ ) { *r1++ = -SYMBOL; *r1++ = *r8; }
3127 }
3128 else {
3129 for ( i = 0; i < -r8[1]; i++ ) {
3130 *r1++ = ARGHEAD+8; *r1++ = 0;
3131 FILLARG(r1);
3132 *r1++ = 8; *r1++ = SYMBOL; *r1++ = 4; *r1++ = *r8;
3133 *r1++ = -1; *r1++ = 1; *r1++ = 1; *r1++ = 3;
3134 }
3135 }
3136 r8 += 2;
3137 }
3138 argout = r1;
3139 }
3140 }
3141/*
3142 #] SYMBOL :
3143 #[ DOTPRODUCT :
3144
3145 Now collect all dotproducts. We can use the space after r1 as storage
3146*/
3147 t = argin2+ARGHEAD;
3148 rnext = t + *t;
3149 r2 = r1;
3150 while ( t < r3 ) {
3151 GETSTOP(t,r6);
3152 t++;
3153 act = 0;
3154 while ( t < r6 ) {
3155 if ( *t == DOTPRODUCT ) {
3156 act = 1;
3157 i = t[1];
3158 NCOPY(r2,t,i)
3159 }
3160 else { t += t[1]; }
3161 }
3162 if ( act == 0 ) {
3163 *r2++ = DOTPRODUCT; *r2++ = 2;
3164 }
3165 t = rnext; rnext = rnext + *rnext;
3166 }
3167 *r2 = 0;
3168 argin3 = argin2;
3169/*
3170 Now we have a list of all dotproducts as a sequence of DOTPRODUCT
3171 subterms. Any dotproduct that is absent in a subterm has power zero.
3172 We now need a list of all minimum powers.
3173 This can be done by subsequent merges.
3174*/
3175 r7 = r1; /* The first object into which we merge. */
3176 r8 = r7 + r7[1]; /* The object that gets merged into r7. */
3177 while ( *r8 ) {
3178 r2 = r8 + r8[1]; /* Next object */
3179 ArgDotproductMerge(r7,r8);
3180 r8 = r2;
3181 }
3182/*
3183 Now we have to divide by the object in r7 and take it apart as factors.
3184 The division can be simple if there are no negative powers.
3185*/
3186 if ( r7[1] > 2 ) {
3187 r8 = r7+2;
3188 r2 = r7 + r7[1];
3189 act = 0;
3190 pow = 0;
3191 while ( r8 < r2 ) {
3192 if ( r8[2] < 0 ) { pow += -r8[2]*(ARGHEAD+9); }
3193 else { pow += r8[2]*(ARGHEAD+9); }
3194 r8 += 3;
3195 }
3196/*
3197 The amount of space we need to move r7 is given in pow
3198 For dotproducts we always need a new location
3199*/
3200 {
3201 argin3 = TermMalloc("TakeArgContent3");
3202/*
3203 We have to multiply the inverse of r7 into argin
3204 The answer should go to argin2.
3205*/
3206 r5 = argin3; *r5++ = 0; *r5++ = 0; FILLARG(r5);
3207 t = argin2+ARGHEAD;
3208 while ( t < r3 ) {
3209 rnext = t + *t;
3210 GETSTOP(t,r6);
3211 r9 = r5;
3212 *r5++ = *t++ + r7[1];
3213 while ( t < r6 ) *r5++ = *t++;
3214 i = r7[1] - 2; r8 = r7+2;
3215 *r5++ = r7[0]; *r5++ = r7[1];
3216 while ( i > 0 ) { *r5++ = *r8++; *r5++ = *r8++; *r5++ = -*r8++; i -= 3; }
3217 while ( t < rnext ) *r5++ = *t++;
3218 Normalize(BHEAD r9);
3219 r5 = r9 + *r9;
3220 }
3221 *r5 = 0;
3222 *argin3 = r5-argin3;
3223/*
3224 We may have to sort the terms in argin3.
3225*/
3226 NewSort(BHEAD0);
3227 t = argin3+ARGHEAD;
3228 while ( *t ) {
3229 StoreTerm(BHEAD t);
3230 t += *t;
3231 }
3232 t = argin3+ARGHEAD;
3233 if ( EndSort(BHEAD t,0) < 0 ) goto Irreg;
3234 while ( *t ) t += *t;
3235 *argin3 = t - argin3;
3236 r3 = t;
3237/*
3238 And now the factors that go to argout.
3239 First we have to move r7 out of the way.
3240*/
3241 r8 = r7+pow; i = r7[1];
3242 while ( --i >= 0 ) r8[i] = r7[i];
3243 r2 += pow;
3244 r8 += 2;
3245 while ( r8 < r2 ) {
3246 for ( i = ABS(r8[2]); i > 0; i-- ) {
3247 *r1++ = ARGHEAD+9; *r1++ = 0; FILLARG(r1);
3248 *r1++ = 9; *r1++ = DOTPRODUCT; *r1++ = 5; *r1++ = *r8;
3249 *r1++ = r8[1]; *r1++ = r8[2] < 0 ? -1: 1;
3250 *r1++ = 1; *r1++ = 1; *r1++ = 3;
3251 }
3252 r8 += 3;
3253 }
3254 argout = r1;
3255 }
3256 }
3257/*
3258 #] DOTPRODUCT :
3259
3260 We have now in argin3 the argument stripped of negative powers and
3261 common factors. The only thing left to deal with is to make the
3262 coefficients integer. For that we have to find the LCM of the denominators
3263 and the GCD of the numerators. And to start with, the sign.
3264 We force the sign of the first term to be positive.
3265*/
3266 t = argin3 + ARGHEAD; pow = 1;
3267 t += *t;
3268 if ( t[-1] < 0 ) {
3269 pow = 0;
3270 t[-1] = -t[-1];
3271 while ( t < r3 ) {
3272 t += *t; t[-1] = -t[-1];
3273 }
3274 }
3275/*
3276 Now the GCD of the numerators and the LCM of the denominators:
3277*/
3278 argfree = TermMalloc("TakeArgContent1");
3279 if ( AN.cmod != 0 ) {
3280 r1 = MakeMod(BHEAD argin3,r1,argfree);
3281 }
3282 else {
3283 r1 = MakeInteger(BHEAD argin3,r1,argfree);
3284 }
3285 if ( pow == 0 ) {
3286 *r1++ = -SNUMBER;
3287 *r1++ = -1;
3288 }
3289 *r1 = 0;
3290/*
3291 Cleanup
3292*/
3293 if ( argin3 != argin2 ) TermFree(argin3,"TakeArgContent3");
3294 if ( argin2 != argin ) TermFree(argin2,"TakeArgContent2");
3295 return(argfree);
3296Irreg:
3297 MesPrint("Irregularity while sorting argument in TakeArgContent");
3298 if ( argin3 != argin2 ) TermFree(argin3,"TakeArgContent3");
3299 if ( argin2 != argin ) TermFree(argin2,"TakeArgContent2");
3300 Terminate(-1);
3301 return(0);
3302}
3303
3304/*
3305 #] TakeArgContent :
3306 #[ MakeInteger :
3307*/
3318WORD *MakeInteger(PHEAD WORD *argin,WORD *argout,WORD *argfree)
3319{
3320 UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMb, *LCMc;
3321 WORD *r, *r1, *r2, *r3, *r4, *r5, *rnext, i, k, j;
3322 WORD kGCD, kLCM, kGCD2, kkLCM, jLCM, jGCD;
3323 GCDbuffer = NumberMalloc("MakeInteger");
3324 GCDbuffer2 = NumberMalloc("MakeInteger");
3325 LCMbuffer = NumberMalloc("MakeInteger");
3326 LCMb = NumberMalloc("MakeInteger");
3327 LCMc = NumberMalloc("MakeInteger");
3328 r4 = argin + *argin;
3329 r = argin + ARGHEAD;
3330/*
3331 First take the first term to load up the LCM and the GCD
3332*/
3333 r2 = r + *r;
3334 j = r2[-1];
3335 r3 = r2 - ABS(j);
3336 k = REDLENG(j);
3337 if ( k < 0 ) k = -k;
3338 while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
3339 for ( kGCD = 0; kGCD < k; kGCD++ ) GCDbuffer[kGCD] = r3[kGCD];
3340 k = REDLENG(j);
3341 if ( k < 0 ) k = -k;
3342 r3 += k;
3343 while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
3344 for ( kLCM = 0; kLCM < k; kLCM++ ) LCMbuffer[kLCM] = r3[kLCM];
3345 r1 = r2;
3346/*
3347 Now go through the rest of the terms in this argument.
3348*/
3349 while ( r1 < r4 ) {
3350 r2 = r1 + *r1;
3351 j = r2[-1];
3352 r3 = r2 - ABS(j);
3353 k = REDLENG(j);
3354 if ( k < 0 ) k = -k;
3355 while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
3356 if ( ( ( GCDbuffer[0] == 1 ) && ( kGCD == 1 ) ) ) {
3357/*
3358 GCD is already 1
3359*/
3360 }
3361 else if ( ( ( k != 1 ) || ( r3[0] != 1 ) ) ) {
3362 if ( GcdLong(BHEAD GCDbuffer,kGCD,(UWORD *)r3,k,GCDbuffer2,&kGCD2) ) {
3363 NumberFree(GCDbuffer,"MakeInteger");
3364 NumberFree(GCDbuffer2,"MakeInteger");
3365 NumberFree(LCMbuffer,"MakeInteger");
3366 NumberFree(LCMb,"MakeInteger"); NumberFree(LCMc,"MakeInteger");
3367 goto MakeIntegerErr;
3368 }
3369 kGCD = kGCD2;
3370 for ( i = 0; i < kGCD; i++ ) GCDbuffer[i] = GCDbuffer2[i];
3371 }
3372 else {
3373 kGCD = 1; GCDbuffer[0] = 1;
3374 }
3375 k = REDLENG(j);
3376 if ( k < 0 ) k = -k;
3377 r3 += k;
3378 while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
3379 if ( ( ( LCMbuffer[0] == 1 ) && ( kLCM == 1 ) ) ) {
3380 for ( kLCM = 0; kLCM < k; kLCM++ )
3381 LCMbuffer[kLCM] = r3[kLCM];
3382 }
3383 else if ( ( k != 1 ) || ( r3[0] != 1 ) ) {
3384 if ( GcdLong(BHEAD LCMbuffer,kLCM,(UWORD *)r3,k,LCMb,&kkLCM) ) {
3385 NumberFree(GCDbuffer,"MakeInteger"); NumberFree(GCDbuffer2,"MakeInteger");
3386 NumberFree(LCMbuffer,"MakeInteger"); NumberFree(LCMb,"MakeInteger"); NumberFree(LCMc,"MakeInteger");
3387 goto MakeIntegerErr;
3388 }
3389 DivLong((UWORD *)r3,k,LCMb,kkLCM,LCMb,&kkLCM,LCMc,&jLCM);
3390 MulLong(LCMbuffer,kLCM,LCMb,kkLCM,LCMc,&jLCM);
3391 for ( kLCM = 0; kLCM < jLCM; kLCM++ )
3392 LCMbuffer[kLCM] = LCMc[kLCM];
3393 }
3394 else {} /* LCM doesn't change */
3395 r1 = r2;
3396 }
3397/*
3398 Now put the factor together: GCD/LCM
3399*/
3400 r3 = (WORD *)(GCDbuffer);
3401 if ( kGCD == kLCM ) {
3402 for ( jGCD = 0; jGCD < kGCD; jGCD++ )
3403 r3[jGCD+kGCD] = LCMbuffer[jGCD];
3404 k = kGCD;
3405 }
3406 else if ( kGCD > kLCM ) {
3407 for ( jGCD = 0; jGCD < kLCM; jGCD++ )
3408 r3[jGCD+kGCD] = LCMbuffer[jGCD];
3409 for ( jGCD = kLCM; jGCD < kGCD; jGCD++ )
3410 r3[jGCD+kGCD] = 0;
3411 k = kGCD;
3412 }
3413 else {
3414 for ( jGCD = kGCD; jGCD < kLCM; jGCD++ )
3415 r3[jGCD] = 0;
3416 for ( jGCD = 0; jGCD < kLCM; jGCD++ )
3417 r3[jGCD+kLCM] = LCMbuffer[jGCD];
3418 k = kLCM;
3419 }
3420 j = 2*k+1;
3421/*
3422 Now we have to write this to argout
3423*/
3424 if ( ( j == 3 ) && ( r3[1] == 1 ) && ( (WORD)(r3[0]) > 0 ) ) {
3425 *argout = -SNUMBER;
3426 argout[1] = r3[0];
3427 r1 = argout+2;
3428 }
3429 else {
3430 r1 = argout;
3431 *r1++ = j+1+ARGHEAD; *r1++ = 0; FILLARG(r1);
3432 *r1++ = j+1; r2 = r3;
3433 for ( i = 0; i < k; i++ ) { *r1++ = *r2++; *r1++ = *r2++; }
3434 *r1++ = j;
3435 }
3436/*
3437 Next we have to take the factor out from the argument.
3438 This cannot be done in location, because the denominator stuff can make
3439 coefficients longer.
3440*/
3441 r2 = argfree + 2; FILLARG(r2)
3442 while ( r < r4 ) {
3443 rnext = r + *r;
3444 j = ABS(rnext[-1]);
3445 r5 = rnext - j;
3446 r3 = r2;
3447 while ( r < r5 ) *r2++ = *r++;
3448 j = (j-1)/2; /* reduced length. Remember, k is the other red length */
3449 if ( DivRat(BHEAD (UWORD *)r5,j,GCDbuffer,k,(UWORD *)r2,&i) ) {
3450 goto MakeIntegerErr;
3451 }
3452 i = 2*i+1;
3453 r2 = r2 + i;
3454 if ( rnext[-1] < 0 ) r2[-1] = -i;
3455 else r2[-1] = i;
3456 *r3 = r2-r3;
3457 r = rnext;
3458 }
3459 *r2 = 0;
3460 argfree[0] = r2-argfree;
3461 argfree[1] = 0;
3462/*
3463 Cleanup
3464*/
3465 NumberFree(LCMc,"MakeInteger");
3466 NumberFree(LCMb,"MakeInteger");
3467 NumberFree(LCMbuffer,"MakeInteger");
3468 NumberFree(GCDbuffer2,"MakeInteger");
3469 NumberFree(GCDbuffer,"MakeInteger");
3470 return(r1);
3471
3472MakeIntegerErr:
3473 MesCall("MakeInteger");
3474 Terminate(-1);
3475 return(0);
3476}
3477
3478/*
3479 #] MakeInteger :
3480 #[ MakeMod :
3481*/
3489WORD *MakeMod(PHEAD WORD *argin,WORD *argout,WORD *argfree)
3490{
3491 WORD *r, *instop, *r1, *m, x, xx, ix, ip;
3492 int i;
3493 r = argin; instop = r + *r; r += ARGHEAD;
3494 x = r[*r-3];
3495 if ( r[*r-1] < 0 ) x += AN.cmod[0];
3496 if ( GetModInverses(x,(WORD)(AN.cmod[0]),&ix,&ip) ) {
3497 Terminate(-1);
3498 }
3499 argout[0] = -SNUMBER;
3500 argout[1] = x;
3501 argout[2] = 0;
3502 r1 = argout+2;
3503/*
3504 Now we have to multiply all coefficients by ix.
3505 This does not make things longer, but we should keep to the conventions
3506 of MakeInteger.
3507*/
3508 m = argfree + ARGHEAD;
3509 while ( r < instop ) {
3510 xx = r[*r-3]; if ( r[*r-1] < 0 ) xx += AN.cmod[0];
3511 xx = (WORD)((((LONG)xx)*ix) % AN.cmod[0]);
3512 if ( xx != 0 ) {
3513 i = *r; NCOPY(m,r,i);
3514 m[-3] = xx; m[-1] = 3;
3515 }
3516 else { r += *r; }
3517 }
3518 *m = 0;
3519 *argfree = m - argfree;
3520 argfree[1] = 0;
3521 argfree += 2; FILLARG(argfree);
3522 return(r1);
3523}
3524
3525/*
3526 #] MakeMod :
3527 #[ SortWeights :
3528*/
3534void SortWeights(LONG *weights,LONG *extraspace,WORD number)
3535{
3536 LONG w, *fill, *from1, *from2;
3537 int n1,n2,i;
3538 if ( number >= 4 ) {
3539 n1 = number/2; n2 = number - n1;
3540 SortWeights(weights,extraspace,n1);
3541 SortWeights(weights+n1,extraspace,n2);
3542/*
3543 We copy the first patch to the extra space. Then we merge
3544 Note that a potential remaining n2 objects are already in place.
3545*/
3546 for ( i = 0; i < n1; i++ ) extraspace[i] = weights[i];
3547 fill = weights; from1 = extraspace; from2 = weights+n1;
3548 while ( n1 > 0 && n2 > 0 ) {
3549 if ( *from1 <= *from2 ) { *fill++ = *from1++; n1--; }
3550 else { *fill++ = *from2++; n2--; }
3551 }
3552 while ( n1 > 0 ) { *fill++ = *from1++; n1--; }
3553 }
3554/*
3555 Special cases
3556*/
3557 else if ( number == 3 ) { /* 6 permutations of which one is trivial */
3558 if ( weights[0] > weights[1] ) {
3559 if ( weights[1] > weights[2] ) {
3560 w = weights[0]; weights[0] = weights[2]; weights[2] = w;
3561 }
3562 else if ( weights[0] > weights[2] ) {
3563 w = weights[0]; weights[0] = weights[1];
3564 weights[1] = weights[2]; weights[2] = w;
3565 }
3566 else {
3567 w = weights[0]; weights[0] = weights[1]; weights[1] = w;
3568 }
3569 }
3570 else if ( weights[0] > weights[2] ) {
3571 w = weights[0]; weights[0] = weights[2];
3572 weights[2] = weights[1]; weights[1] = w;
3573 }
3574 else if ( weights[1] > weights[2] ) {
3575 w = weights[1]; weights[1] = weights[2]; weights[2] = w;
3576 }
3577 }
3578 else if ( number == 2 ) {
3579 if ( weights[0] > weights[1] ) {
3580 w = weights[0]; weights[0] = weights[1]; weights[1] = w;
3581 }
3582 }
3583 return;
3584}
3585
3586/*
3587 #] SortWeights :
3588*/
TakeArgContent
WORD * TakeArgContent(PHEAD WORD *argin, WORD *argout)
Definition argument.c:2772
MakeMod
WORD * MakeMod(PHEAD WORD *argin, WORD *argout, WORD *argfree)
Definition argument.c:3489
CleanupArgCache
int CleanupArgCache(PHEAD WORD bufnum)
Definition argument.c:2578
MakeInteger
WORD * MakeInteger(PHEAD WORD *argin, WORD *argout, WORD *argfree)
Definition argument.c:3318
InsertArg
int InsertArg(PHEAD WORD *argin, WORD *argout, int par)
Definition argument.c:2543
SortWeights
void SortWeights(LONG *weights, LONG *extraspace, WORD number)
Definition argument.c:3534
FindArg
WORD FindArg(PHEAD WORD *a)
Definition argument.c:2519
AddRHS
WORD * AddRHS(int num, int type)
Definition comtool.c:214
DoubleCbuffer
WORD * DoubleCbuffer(int num, WORD *w, int par)
Definition comtool.c:143
AddNtoC
int AddNtoC(int bufnum, int n, WORD *array, int par)
Definition comtool.c:317
LocalConvertToPoly
int LocalConvertToPoly(PHEAD WORD *, WORD *, WORD, WORD)
Definition notation.c:510
EndSort
LONG EndSort(PHEAD WORD *, int)
Definition sort.c:454
ZeroFillRange
#define ZeroFillRange(w, begin, end)
Definition declare.h:182
Generator
int Generator(PHEAD WORD *, WORD)
Definition proces.c:3249
LowerSortLevel
void LowerSortLevel(void)
Definition sort.c:4661
StoreTerm
int StoreTerm(PHEAD WORD *)
Definition sort.c:4244
NewSort
int NewSort(PHEAD0)
Definition sort.c:359
poly_factorize_argument
int poly_factorize_argument(PHEAD WORD *, WORD *)
Definition polywrap.cc:1112
GetModInverses
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition reken.c:1477
CbUf::boomlijst
COMPTREE * boomlijst
Definition structs.h:980
CbUf::rhs
WORD ** rhs
Definition structs.h:975
CbUf::lhs
WORD ** lhs
Definition structs.h:974
CbUf::Buffer
WORD * Buffer
Definition structs.h:971
CbUf::Pointer
WORD * Pointer
Definition structs.h:973
tree::usage
int usage
Definition structs.h:294
Generated by doxygen 1.9.8