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