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