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