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