FORM v5.0.2-2-gfb12f2e
normal.c
Go to the documentation of this file.
1
11/* #[ License : */
12/*
13 * Copyright (C) 1984-2026 J.A.M. Vermaseren
14 * When using this file you are requested to refer to the publication
15 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
16 * This is considered a matter of courtesy as the development was paid
17 * for by FOM the Dutch physics granting agency and we would like to
18 * be able to track its scientific use to convince FOM of its value
19 * for the community.
20 *
21 * This file is part of FORM.
22 *
23 * FORM is free software: you can redistribute it and/or modify it under the
24 * terms of the GNU General Public License as published by the Free Software
25 * Foundation, either version 3 of the License, or (at your option) any later
26 * version.
27 *
28 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
29 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
30 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
31 * details.
32 *
33 * You should have received a copy of the GNU General Public License along
34 * with FORM. If not, see <http://www.gnu.org/licenses/>.
35 */
36/* #] License : */
37/*
38 #[ Includes : normal.c
39*/
40
41#include "form3.h"
42#ifdef WITHFLOAT
43#include <gmp.h>
44
45int PackFloat(WORD *,mpf_t);
46int UnpackFloat(mpf_t, WORD *);
47void RatToFloat(mpf_t result, UWORD *formrat, int ratsize);
48#endif
49/*
50 #] Includes :
51 #[ Normalize :
52 #[ CompareFunctions :
53*/
54
55int CompareFunctions(WORD *fleft,WORD *fright)
56{
57 WORD k, kk;
58 if ( AC.properorderflag ) {
59 if ( ( *fleft >= (FUNCTION+WILDOFFSET)
60 && functions[*fleft-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
61 || ( *fleft >= FUNCTION && *fleft < (FUNCTION + WILDOFFSET)
62 && functions[*fleft-FUNCTION].spec >= TENSORFUNCTION ) ) {}
63 else {
64 WORD *s1, *s2, *ss1, *ss2;
65 s1 = fleft+FUNHEAD; s2 = fright+FUNHEAD;
66 ss1 = fleft + fleft[1]; ss2 = fright + fright[1];
67 while ( s1 < ss1 && s2 < ss2 ) {
68 k = CompArg(s1,s2);
69 if ( k > 0 ) return(1);
70 if ( k < 0 ) return(0);
71 NEXTARG(s1)
72 NEXTARG(s2)
73 }
74 if ( s1 < ss1 ) return(1);
75 return(0);
76 }
77 k = fleft[1] - FUNHEAD;
78 kk = fright[1] - FUNHEAD;
79 fleft += FUNHEAD;
80 fright += FUNHEAD;
81 while ( k > 0 && kk > 0 ) {
82 if ( *fleft < *fright ) return(0);
83 else if ( *fleft++ > *fright++ ) return(1);
84 k--; kk--;
85 }
86 if ( k > 0 ) return(1);
87 return(0);
88 }
89 else {
90 k = fleft[1] - FUNHEAD;
91 kk = fright[1] - FUNHEAD;
92 fleft += FUNHEAD;
93 fright += FUNHEAD;
94 while ( k > 0 && kk > 0 ) {
95 if ( *fleft < *fright ) return(0);
96 else if ( *fleft++ > *fright++ ) return(1);
97 k--; kk--;
98 }
99 if ( k > 0 ) return(1);
100 return(0);
101 }
102}
103
104/*
105 #] CompareFunctions :
106 #[ Commute :
107
108 This function gets two adjacent function pointers and decides
109 whether these two functions should be exchanged to obtain a
110 natural ordering.
111
112 Currently there is only an ordering of gamma matrices belonging
113 to different spin lines.
114
115 Note that we skip for now the cases of (F)^(3/2) or 1/F and a few more
116 of such funny functions.
117*/
118
119int Commute(WORD *fleft, WORD *fright)
120{
121 WORD fun1, fun2;
122 if ( *fleft == DOLLAREXPRESSION || *fright == DOLLAREXPRESSION ) return(0);
123 fun1 = ABS(*fleft); fun2 = ABS(*fright);
124 if ( *fleft >= GAMMA && *fleft <= GAMMASEVEN
125 && *fright >= GAMMA && *fright <= GAMMASEVEN ) {
126 if ( fleft[FUNHEAD] < AM.OffsetIndex && fleft[FUNHEAD] > fright[FUNHEAD] )
127 return(1);
128 return(0);
129 }
130 if ( fun1 >= WILDOFFSET ) fun1 -= WILDOFFSET;
131 if ( fun2 >= WILDOFFSET ) fun2 -= WILDOFFSET;
132 if ( ( ( functions[fun1-FUNCTION].flags & COULDCOMMUTE ) == 0 )
133 || ( ( functions[fun2-FUNCTION].flags & COULDCOMMUTE ) == 0 ) ) return(0);
134/*
135 if other conditions will come here, keep in mind that if *fleft < 0
136 or *fright < 0 they are arguments in the exponent function as in f^(3/2)
137*/
138 if ( AC.CommuteInSet == 0 ) return(0);
139/*
140 The code for CompareFunctions can be stolen from the commuting case.
141
142 We need the syntax:
143 Commute Fun1,Fun2,...,Fun`n';
144 For this Fun1,...,Fun`n' need to be noncommuting functions.
145 These functions will commute with all members of the group.
146 In the AC.paircommute buffer the representation is
147 `n'+1,element1,...,element`n',`m'+1,element1,...,element`m',0
148 A function can belong to more than one group.
149 If a function commutes with itself, it is most efficient to make a separate
150 group of two elements for it as in
151 Commute T,T; -> 3,T,T
152*/
153 if ( fun1 >= fun2 ) {
154 WORD *group = AC.CommuteInSet, *g1, *g2, *g3;
155 while ( *group > 0 ) {
156 g3 = group + *group;
157 g1 = group+1;
158 while ( g1 < g3 ) {
159 if ( *g1 == fun1 || ( fun1 <= GAMMASEVEN && fun1 >= GAMMA
160 && *g1 <= GAMMASEVEN && *g1 >= GAMMA ) ) {
161 g2 = group+1;
162 while ( g2 < g3 ) {
163 if ( g1 != g2 && ( *g2 == fun2 ||
164 ( fun2 <= GAMMASEVEN && fun2 >= GAMMA
165 && *g2 <= GAMMASEVEN && *g2 >= GAMMA ) ) ) {
166 if ( fun1 != fun2 ) return(1);
167 if ( *fleft < 0 ) return(0);
168 if ( *fright < 0 ) return(1);
169 return(CompareFunctions(fleft,fright));
170 }
171 g2++;
172 }
173 break;
174 }
175 g1++;
176 }
177 group = g3;
178 }
179 }
180 return(0);
181}
182
183/*
184 #] Commute :
185 #[ Normalize :
186
187 This is the big normalization routine. It has a great need
188 to be economical.
189 The limit on the number of objects coming in is given by NORMSIZE.
190
191*/
192
193int Normalize(PHEAD WORD *term)
194{
195/*
196 #[ Declarations :
197*/
198 GETBIDENTITY
199 WORD *t, *m, *r, i, j, k, l, nsym, *ss, *tt, *u;
200 WORD shortnum, stype;
201 WORD *stop, *to = 0, *from = 0;
202
203 WORD *ppsym, *ppvec, *ppdot, *ppdel;
204 WORD nvec, ndot, ndel, nind, neps, nden, ncom, nnco, ncon;
205
206 AT.NormDepth++;
207#ifdef DEBUGGING
208 if ( AT.NormDepth > 2 ) {
209 // We don't expect this to happen in the current codebase.
210 Terminate(-1);
211 }
212#endif
213 if ( AT.NormDepth > AT.NormDataSize ) {
214 NORMDATA **top = AT.NormData + AT.NormDataSize;
215 DoubleBuffer((void **)&(AT.NormData), (void **)&(top),
216 sizeof(*AT.NormData), "double NormData pointers");
217 AT.NormDataSize *= 2;
218 for ( LONG i = AT.NormDepth-1; i < AT.NormDataSize; i++ ) {
219 AT.NormData[i] = NULL;
220 }
221 }
222 if ( AT.NormData[AT.NormDepth-1] == NULL ) {
223 AT.NormData[AT.NormDepth-1] = AllocNormData();
224 }
225
226 WORD *psym = AT.NormData[AT.NormDepth-1]->psym;
227 WORD *pvec = AT.NormData[AT.NormDepth-1]->pvec;
228 WORD *pdot = AT.NormData[AT.NormDepth-1]->pdot;
229 WORD *pdel = AT.NormData[AT.NormDepth-1]->pdel;
230 WORD *pind = AT.NormData[AT.NormDepth-1]->pind;
231 WORD **peps = AT.NormData[AT.NormDepth-1]->peps;
232 WORD **pden = AT.NormData[AT.NormDepth-1]->pden;
233 WORD **pcom = AT.NormData[AT.NormDepth-1]->pcom;
234 WORD **pnco = AT.NormData[AT.NormDepth-1]->pnco;
235 WORD **pcon = AT.NormData[AT.NormDepth-1]->pcon;
236
237 WORD *n_coef, ncoef; /* Accumulator for the coefficient */
238 WORD *n_llnum, *lnum, nnum;
239 WORD *termout, oldtoprhs = 0, subtype;
240 WORD ReplaceType, ReplaceVeto = 0, didcontr;
241 int regval = 0;
242 WORD *ReplaceSub;
243 WORD *fillsetexp;
244 CBUF *C = cbuf+AT.ebufnum;
245 WORD *ANsc = 0, *ANsm = 0, *ANsr = 0, PolyFunMode;
246#ifdef WITHFLOAT
247 WORD withfloat = 0;
248 WORD *firstfloat = 0;
249#endif
250 LONG oldcpointer = 0, x;
251 n_coef = TermMalloc("NormCoef");
252 n_llnum = TermMalloc("n_llnum");
253 lnum = n_llnum+1;
254/*
255 int termflag;
256*/
257/*
258 #] Declarations :
259 #[ Setup :
260PrintTerm(term,"Normalize");
261*/
262
263 if ( AT.SS == AT.S0 ) {
264 if ( *term > AT.SS->verbMaxTermSize ) {
265 AT.SS->verbMaxTermSize = *term;
266 }
267 }
268
269Restart:
270 didcontr = 0;
271 ReplaceType = -1;
272 t = term;
273 if ( !*t ) {
274 AT.NormDepth--;
275 TermFree(n_coef,"NormCoef");
276 TermFree(n_llnum,"n_llnum");
277 return(regval);
278 }
279 r = t + *t;
280 ncoef = r[-1];
281 i = ABS(ncoef);
282 r -= i;
283 m = r;
284 t = n_coef;
285 NCOPY(t,r,i);
286 termout = AT.WorkPointer;
287 AT.WorkPointer = (WORD *)(((UBYTE *)(AT.WorkPointer)) + AM.MaxTer);
288 fillsetexp = termout+1;
289 AN.PolyNormFlag = 0; PolyFunMode = AN.PolyFunTodo;
290/*
291 termflag = 0;
292*/
293/*
294 #] Setup :
295 #[ First scan :
296*/
297 nsym = nvec = ndot = ndel = neps = nden =
298 nind = ncom = nnco = ncon = 0;
299 ppsym = psym;
300 ppvec = pvec;
301 ppdot = pdot;
302 ppdel = pdel;
303 t = term + 1;
304conscan:;
305 if ( t < m ) do {
306 r = t + t[1];
307 switch ( *t ) {
308 case SYMBOL :
309 t += 2;
310 from = m;
311 do {
312 if ( t[1] == 0 ) {
313/* if ( *t == 0 || *t == MAXPOWER ) goto NormZZ; */
314 t += 2;
315 goto NextSymbol;
316 }
317 if ( *t <= DENOMINATORSYMBOL && *t >= COEFFSYMBOL ) {
318/*
319 if ( AN.NoScrat2 == 0 ) {
320 AN.NoScrat2 = (UWORD *)Malloc1((AM.MaxTal+2)*sizeof(UWORD),"Normalize");
321 }
322*/
323 if ( AN.cTerm ) m = AN.cTerm;
324 else m = term;
325 m += *m;
326 ncoef = REDLENG(ncoef);
327 if ( *t == COEFFSYMBOL ) {
328 i = t[1];
329 nnum = REDLENG(m[-1]);
330 m -= ABS(m[-1]);
331 if ( i > 0 ) {
332 while ( i > 0 ) {
333 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,nnum,
334 (UWORD *)n_coef,&ncoef) ) goto FromNorm;
335 i--;
336 }
337 }
338 else if ( i < 0 ) {
339 while ( i < 0 ) {
340 if ( DivRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,nnum,
341 (UWORD *)n_coef,&ncoef) ) goto FromNorm;
342 i++;
343 }
344 }
345 }
346 else {
347 i = m[-1];
348 nnum = (ABS(i)-1)/2;
349 if ( *t == NUMERATORSYMBOL ) { m -= nnum + 1; }
350 else { m--; }
351 while ( *m == 0 && nnum > 1 ) { m--; nnum--; }
352 m -= nnum;
353 if ( i < 0 && *t == NUMERATORSYMBOL ) nnum = -nnum;
354 i = t[1];
355 if ( i > 0 ) {
356 while ( i > 0 ) {
357 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)m,nnum) )
358 goto FromNorm;
359 i--;
360 }
361 }
362 else if ( i < 0 ) {
363 while ( i < 0 ) {
364 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)m,nnum) )
365 goto FromNorm;
366 i++;
367 }
368 }
369 }
370 ncoef = INCLENG(ncoef);
371 t += 2;
372 goto NextSymbol;
373 }
374 else if ( *t == DIMENSIONSYMBOL ) {
375 if ( AN.cTerm ) m = AN.cTerm;
376 else m = term;
377 k = DimensionTerm(m);
378 if ( k == 0 ) goto NormZero;
379 if ( k == MAXPOSITIVE ) {
380 MLOCK(ErrorMessageLock);
381 MesPrint("Dimension_ is undefined in term %t");
382 MUNLOCK(ErrorMessageLock);
383 goto NormMin;
384 }
385 if ( k == -MAXPOSITIVE ) {
386 MLOCK(ErrorMessageLock);
387 MesPrint("Dimension_ out of range in term %t");
388 MUNLOCK(ErrorMessageLock);
389 goto NormMin;
390 }
391 if ( k > 0 ) { *((UWORD *)lnum) = k; nnum = 1; }
392 else { *((UWORD *)lnum) = -k; nnum = -1; }
393 ncoef = REDLENG(ncoef);
394 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) goto FromNorm;
395 ncoef = INCLENG(ncoef);
396 t += 2;
397 goto NextSymbol;
398 }
399 if ( ( *t >= MAXPOWER && *t < 2*MAXPOWER )
400 || ( *t < -MAXPOWER && *t > -2*MAXPOWER ) ) {
401/*
402 #[ TO SNUMBER :
403*/
404 if ( *t < 0 ) {
405 *t += MAXPOWER;
406 *t = -*t;
407 if ( t[1] & 1 ) ncoef = -ncoef;
408 }
409 else if ( *t == MAXPOWER ) {
410 if ( t[1] > 0 ) goto NormZero;
411 goto NormInf;
412 }
413 else {
414 *t -= MAXPOWER;
415 }
416 lnum[0] = *t;
417 nnum = 1;
418 if ( t[1] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[1]))) )
419 goto FromNorm;
420 ncoef = REDLENG(ncoef);
421 if ( t[1] < 0 ) {
422 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
423 goto FromNorm;
424 }
425 else if ( t[1] > 0 ) {
426 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
427 goto FromNorm;
428 }
429 ncoef = INCLENG(ncoef);
430/*
431 #] TO SNUMBER :
432*/
433 t += 2;
434 goto NextSymbol;
435 }
436 if ( ( *t <= NumSymbols && *t > -MAXPOWER )
437 && ( symbols[*t].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
438 if ( t[1] <= 2*MAXPOWER && t[1] >= -2*MAXPOWER ) {
439 t[1] %= symbols[*t].maxpower;
440 if ( t[1] < 0 ) t[1] += symbols[*t].maxpower;
441 if ( ( symbols[*t].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
442 if ( ( ( symbols[*t].maxpower & 1 ) == 0 ) &&
443 ( t[1] >= symbols[*t].maxpower/2 ) ) {
444 t[1] -= symbols[*t].maxpower/2; ncoef = -ncoef;
445 }
446 }
447 if ( t[1] == 0 ) { t += 2; goto NextSymbol; }
448 }
449 }
450 i = nsym;
451 m = ppsym;
452 if ( i > 0 ) do {
453 m -= 2;
454 if ( *t == *m ) {
455 t++; m++;
456 if ( *t > 2*MAXPOWER || *t < -2*MAXPOWER
457 || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
458 MLOCK(ErrorMessageLock);
459 MesPrint("Illegal wildcard power combination.");
460 MUNLOCK(ErrorMessageLock);
461 goto NormMin;
462 }
463 *m += *t;
464 if ( ( t[-1] <= NumSymbols && t[-1] > -MAXPOWER )
465 && ( symbols[t[-1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
466 *m %= symbols[t[-1]].maxpower;
467 if ( *m < 0 ) *m += symbols[t[-1]].maxpower;
468 if ( ( symbols[t[-1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
469 if ( ( ( symbols[t[-1]].maxpower & 1 ) == 0 ) &&
470 ( *m >= symbols[t[-1]].maxpower/2 ) ) {
471 *m -= symbols[t[-1]].maxpower/2; ncoef = -ncoef;
472 }
473 }
474 }
475 if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
476 MLOCK(ErrorMessageLock);
477 MesPrint("Power overflow during normalization");
478 MUNLOCK(ErrorMessageLock);
479 goto NormMin;
480 }
481 if ( !*m ) {
482 m--;
483 while ( i < nsym )
484 { *m = m[2]; m++; *m = m[2]; m++; i++; }
485 ppsym -= 2;
486 nsym--;
487 }
488 t++;
489 goto NextSymbol;
490 }
491 } while ( *t < *m && --i > 0 );
492 m = ppsym;
493 while ( i < nsym )
494 { m--; m[2] = *m; m--; m[2] = *m; i++; }
495 *m++ = *t++;
496 *m = *t++;
497 ppsym += 2;
498 nsym++;
499NextSymbol:;
500 } while ( t < r );
501 m = from;
502 break;
503 case VECTOR :
504 t += 2;
505 do {
506 if ( t[0] == AM.vectorzero ) goto NormZero;
507 if ( t[1] == FUNNYVEC ) {
508 pind[nind++] = *t;
509 t += 2;
510 }
511 else if ( t[1] < 0 ) {
512 if ( *t == NOINDEX && t[1] == NOINDEX ) t += 2;
513 else {
514 if ( t[1] == AM.vectorzero ) goto NormZero;
515 *ppdot++ = *t++; *ppdot++ = *t++; *ppdot++ = 1; ndot++;
516 }
517 }
518 else { *ppvec++ = *t++; *ppvec++ = *t++; nvec += 2; }
519 } while ( t < r );
520 break;
521 case DOTPRODUCT :
522 t += 2;
523 do {
524 if ( t[2] == 0 ) t += 3;
525 else if ( ndot > 0 && t[0] == ppdot[-3]
526 && t[1] == ppdot[-2] ) {
527 ppdot[-1] += t[2];
528 t += 3;
529 if ( ppdot[-1] == 0 ) { ppdot -= 3; ndot--; }
530 }
531 else if ( t[0] == AM.vectorzero || t[1] == AM.vectorzero ) {
532 if ( t[2] > 0 ) goto NormZero;
533 goto NormInf;
534 }
535 else {
536 *ppdot++ = *t++; *ppdot++ = *t++;
537 *ppdot++ = *t++; ndot++;
538 }
539 } while ( t < r );
540 break;
541 case HAAKJE :
542 break;
543 case SETSET:
544 if ( WildFill(BHEAD termout,term,AT.dummysubexp) < 0 ) goto FromNorm;
545 i = *termout;
546 t = termout; m = term;
547 NCOPY(m,t,i);
548 goto Restart;
549 case DOLLAREXPRESSION :
550/*
551 We have DOLLAREXPRESSION,4,number,power
552 Replace by SUBEXPRESSION and exit elegantly to let
553 TestSub pick it up. Of course look for special cases first.
554 Note that we have a special compiler buffer for the values.
555*/
556 if ( AR.Eside != LHSIDE ) {
557 DOLLARS d = Dollars + t[2];
558#ifdef WITHPTHREADS
559 int nummodopt, ptype = -1;
560 if ( AS.MultiThreaded ) {
561 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
562 if ( t[2] == ModOptdollars[nummodopt].number ) break;
563 }
564 if ( nummodopt < NumModOptdollars ) {
565 ptype = ModOptdollars[nummodopt].type;
566 if ( DollarLocalCopy(ptype) ) {
567 d = ModOptdollars[nummodopt].dstruct+AT.identity;
568 }
569 else {
570 LOCK(d->pthreadslock);
571 }
572 }
573 }
574#endif
575 if ( d->type == DOLZERO ) {
576#ifdef WITHPTHREADS
577 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
578#endif
579 if ( t[3] == 0 ) goto NormZZ;
580 if ( t[3] < 0 ) goto NormInf;
581 goto NormZero;
582 }
583 else if ( d->type == DOLNUMBER ) {
584 nnum = d->where[0];
585 if ( nnum > 0 ) {
586 nnum = d->where[nnum-1];
587 if ( nnum < 0 ) { ncoef = -ncoef; nnum = -nnum; }
588 nnum = (nnum-1)/2;
589 for ( i = 1; i <= nnum; i++ ) lnum[i-1] = d->where[i];
590 }
591 if ( nnum == 0 || ( nnum == 1 && lnum[0] == 0 ) ) {
592#ifdef WITHPTHREADS
593 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
594#endif
595 if ( t[3] < 0 ) goto NormInf;
596 else if ( t[3] == 0 ) goto NormZZ;
597 goto NormZero;
598 }
599 if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) ) goto FromNorm;
600 ncoef = REDLENG(ncoef);
601 if ( t[3] < 0 ) {
602 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) {
603#ifdef WITHPTHREADS
604 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
605#endif
606 goto FromNorm;
607 }
608 }
609 else if ( t[3] > 0 ) {
610 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) {
611#ifdef WITHPTHREADS
612 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
613#endif
614 goto FromNorm;
615 }
616 }
617 ncoef = INCLENG(ncoef);
618 }
619 else if ( d->type == DOLINDEX ) {
620 if ( d->index == 0 ) {
621#ifdef WITHPTHREADS
622 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
623#endif
624 goto NormZero;
625 }
626 if ( d->index != NOINDEX ) pind[nind++] = d->index;
627 }
628 else if ( d->type == DOLTERMS ) {
629 if ( t[3] >= MAXPOWER || t[3] <= -MAXPOWER ) {
630 if ( d->where[0] == 0 ) goto NormZero;
631 if ( d->where[d->where[0]] != 0 ) {
632IllDollarExp:
633 MLOCK(ErrorMessageLock);
634 MesPrint("!!!Illegal $ expansion with wildcard power!!!");
635 MUNLOCK(ErrorMessageLock);
636 goto FromNorm;
637 }
638/*
639 At this point we should only admit symbols and dotproducts
640 We expand the dollar directly and do not send it back.
641*/
642 { WORD *td, *tdstop, dj;
643 td = d->where+1;
644 tdstop = d->where+d->where[0];
645 if ( tdstop[-1] != 3 || tdstop[-2] != 1
646 || tdstop[-3] != 1 ) goto IllDollarExp;
647 tdstop -= 3;
648 if ( td >= tdstop ) goto IllDollarExp;
649 while ( td < tdstop ) {
650 if ( *td == SYMBOL ) {
651 for ( dj = 2; dj < td[1]; dj += 2 ) {
652 if ( td[dj+1] == 1 ) {
653 *ppsym++ = td[dj];
654 *ppsym++ = t[3];
655 nsym++;
656 }
657 else if ( td[dj+1] == -1 ) {
658 *ppsym++ = td[dj];
659 *ppsym++ = -t[3];
660 nsym++;
661 }
662 else goto IllDollarExp;
663 }
664 }
665 else if ( *td == DOTPRODUCT ) {
666 for ( dj = 2; dj < td[1]; dj += 3 ) {
667 if ( td[dj+2] == 1 ) {
668 *ppdot++ = td[dj];
669 *ppdot++ = td[dj+1];
670 *ppdot++ = t[3];
671 ndot++;
672 }
673 else if ( td[dj+2] == -1 ) {
674 *ppdot++ = td[dj];
675 *ppdot++ = td[dj+1];
676 *ppdot++ = -t[3];
677 ndot++;
678 }
679 else goto IllDollarExp;
680 }
681 }
682 else goto IllDollarExp;
683 td += td[1];
684 }
685 regval = 2;
686 break;
687 }
688 }
689 t[0] = SUBEXPRESSION;
690 t[4] = AM.dbufnum;
691 if ( t[3] == 0 ) {
692#ifdef WITHPTHREADS
693 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
694#endif
695 break;
696 }
697 regval = 2;
698 t = r;
699 while ( t < m ) {
700 if ( *t == DOLLAREXPRESSION ) {
701#ifdef WITHPTHREADS
702 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
703#endif
704 d = Dollars + t[2];
705#ifdef WITHPTHREADS
706 if ( AS.MultiThreaded ) {
707 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
708 if ( t[2] == ModOptdollars[nummodopt].number ) break;
709 }
710 if ( nummodopt < NumModOptdollars ) {
711 ptype = ModOptdollars[nummodopt].type;
712 if ( DollarLocalCopy(ptype) ) {
713 d = ModOptdollars[nummodopt].dstruct+AT.identity;
714 }
715 else {
716 LOCK(d->pthreadslock);
717 }
718 }
719 }
720#endif
721 if ( d->type == DOLTERMS ) {
722 t[0] = SUBEXPRESSION;
723 t[4] = AM.dbufnum;
724 }
725 }
726 t += t[1];
727 }
728#ifdef WITHPTHREADS
729 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
730#endif
731 goto RegEnd;
732 }
733 else {
734#ifdef WITHPTHREADS
735 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
736#endif
737 MLOCK(ErrorMessageLock);
738 MesPrint("!!!This $ variation has not been implemented yet!!!");
739 MUNLOCK(ErrorMessageLock);
740 goto NormMin;
741 }
742#ifdef WITHPTHREADS
743 if ( ptype > 0 && ! DollarLocalCopy(ptype) ) { UNLOCK(d->pthreadslock); }
744#endif
745 }
746 else {
747 pnco[nnco++] = t;
748/*
749 The next statement should be safe as the value is used
750 only by the compiler (ie the master).
751*/
752 AC.lhdollarflag = 1;
753 }
754 break;
755 case DELTA :
756 t += 2;
757 do {
758 if ( *t < 0 ) {
759 if ( *t == SUMMEDIND ) {
760 if ( t[1] < -NMIN4SHIFT ) {
761 k = -t[1]-NMIN4SHIFT;
762 k = ExtraSymbol(k,1,nsym,ppsym,&ncoef);
763 nsym += k;
764 ppsym += (k * 2);
765 }
766 else if ( t[1] == 0 ) goto NormZero;
767 else {
768 if ( t[1] < 0 ) {
769 lnum[0] = -t[1];
770 nnum = -1;
771 }
772 else {
773 lnum[0] = t[1];
774 nnum = 1;
775 }
776 ncoef = REDLENG(ncoef);
777 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
778 goto FromNorm;
779 ncoef = INCLENG(ncoef);
780 }
781 t += 2;
782 }
783 else if ( *t == NOINDEX && t[1] == NOINDEX ) t += 2;
784 else if ( *t == EMPTYINDEX && t[1] == EMPTYINDEX ) {
785 *ppdel++ = *t++; *ppdel++ = *t++; ndel += 2;
786 }
787 else
788 if ( t[1] < 0 ) {
789 *ppdot++ = *t++; *ppdot++ = *t++; *ppdot++ = 1; ndot++;
790 }
791 else {
792 *ppvec++ = *t++; *ppvec++ = *t++; nvec += 2;
793 }
794 }
795 else {
796 if ( t[1] < 0 ) {
797 *ppvec++ = t[1]; *ppvec++ = *t; t+=2; nvec += 2;
798 }
799 else { *ppdel++ = *t++; *ppdel++ = *t++; ndel += 2; }
800 }
801 } while ( t < r );
802 break;
803 case FACTORIAL :
804/*
805 (FACTORIAL,FUNHEAD+2,..,-SNUMBER,number)
806*/
807 if ( t[FUNHEAD] == -SNUMBER && t[1] == FUNHEAD+2
808 && t[FUNHEAD+1] >= 0 ) {
809 if ( Factorial(BHEAD t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
810 goto FromNorm;
811MulIn:
812 ncoef = REDLENG(ncoef);
813 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) goto FromNorm;
814 ncoef = INCLENG(ncoef);
815 }
816 else pcom[ncom++] = t;
817 break;
818 case BERNOULLIFUNCTION :
819/*
820 (BERNOULLIFUNCTION,FUNHEAD+2,..,-SNUMBER,number)
821*/
822 if ( ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] >= 0 )
823 && ( t[1] == FUNHEAD+2 || ( t[1] == FUNHEAD+4 &&
824 t[FUNHEAD+2] == -SNUMBER && ABS(t[FUNHEAD+3]) == 1 ) ) ) {
825 WORD inum, mnum;
826 if ( Bernoulli(t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
827 goto FromNorm;
828 if ( nnum == 0 ) goto NormZero;
829 inum = nnum; if ( inum < 0 ) inum = -inum;
830 inum--; inum /= 2;
831 mnum = inum;
832 while ( lnum[mnum-1] == 0 ) mnum--;
833 if ( nnum < 0 ) mnum = -mnum;
834 ncoef = REDLENG(ncoef);
835 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,mnum) ) goto FromNorm;
836 mnum = inum;
837 while ( lnum[inum+mnum-1] == 0 ) mnum--;
838 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(lnum+inum),mnum) ) goto FromNorm;
839 ncoef = INCLENG(ncoef);
840 if ( t[1] == FUNHEAD+4 && t[FUNHEAD+1] == 1
841 && t[FUNHEAD+3] == -1 ) ncoef = -ncoef;
842 }
843 else pcom[ncom++] = t;
844 break;
845 case NUMARGSFUN:
846/*
847 Numerical function giving the number of arguments.
848*/
849 k = 0;
850 t += FUNHEAD;
851 while ( t < r ) {
852 k++;
853 NEXTARG(t);
854 }
855 if ( k == 0 ) goto NormZero;
856 *((UWORD *)lnum) = k;
857 nnum = 1;
858 goto MulIn;
859 case NUMFACTORS:
860/*
861 Numerical function giving the number of factors in an expression.
862*/
863 t += FUNHEAD;
864 if ( *t == -EXPRESSION ) {
865 k = AS.OldNumFactors[t[1]];
866 }
867 else if ( *t == -DOLLAREXPRESSION ) {
868 k = Dollars[t[1]].nfactors;
869 }
870 else {
871 pcom[ncom++] = t;
872 break;
873 }
874 if ( k == 0 ) goto NormZero;
875 *((UWORD *)lnum) = k;
876 nnum = 1;
877 goto MulIn;
878 case NUMTERMSFUN:
879/*
880 Numerical function giving the number of terms in the single argument.
881*/
882 if ( t[FUNHEAD] < 0 ) {
883 if ( t[FUNHEAD] <= -FUNCTION && t[1] == FUNHEAD+1 ) break;
884 if ( t[FUNHEAD] > -FUNCTION && t[1] == FUNHEAD+2 ) {
885 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] == 0 ) goto NormZero;
886 break;
887 }
888 pcom[ncom++] = t;
889 break;
890 }
891 if ( t[FUNHEAD] > 0 && t[FUNHEAD] == t[1]-FUNHEAD ) {
892 k = 0;
893 t += FUNHEAD+ARGHEAD;
894 while ( t < r ) {
895 k++;
896 t += *t;
897 }
898 if ( k == 0 ) goto NormZero;
899 *((UWORD *)lnum) = k;
900 nnum = 1;
901 goto MulIn;
902 }
903 else pcom[ncom++] = t;
904 break;
905 case COUNTFUNCTION:
906 if ( AN.cTerm ) {
907 k = CountFun(AN.cTerm,t);
908 }
909 else {
910 k = CountFun(term,t);
911 }
912 if ( k == 0 ) goto NormZero;
913 if ( k > 0 ) { *((UWORD *)lnum) = k; nnum = 1; }
914 else { *((UWORD *)lnum) = -k; nnum = -1; }
915 goto MulIn;
916 break;
917 case MAKERATIONAL:
918 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+2] == -SNUMBER
919 && t[1] == FUNHEAD+4 && t[FUNHEAD+3] > 1 ) {
920 WORD x1[2], sgn;
921 if ( t[FUNHEAD+1] == 0 ) goto NormZero;
922 if ( t[FUNHEAD+1] < 0 ) { t[FUNHEAD+1] = -t[FUNHEAD+1]; sgn = -1; }
923 else sgn = 1;
924 if ( MakeRational(t[FUNHEAD+1],t[FUNHEAD+3],x1,x1+1) ) {
925 static int warnflag = 1;
926 if ( warnflag ) {
927 MesPrint("%w Warning: fraction could not be reconstructed in MakeRational_");
928 warnflag = 0;
929 }
930 x1[0] = t[FUNHEAD+1]; x1[1] = 1;
931 }
932 if ( sgn < 0 ) { t[FUNHEAD+1] = -t[FUNHEAD+1]; x1[0] = -x1[0]; }
933 if ( x1[0] < 0 ) { sgn = -1; x1[0] = -x1[0]; }
934 else sgn = 1;
935 ncoef = REDLENG(ncoef);
936 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)x1,sgn,
937 (UWORD *)n_coef,&ncoef) ) goto FromNorm;
938 ncoef = INCLENG(ncoef);
939 }
940 else {
941 WORD narg = 0, *tt, *ttstop, *arg1 = 0, *arg2 = 0;
942 UWORD *x1, *x2, *xx;
943 WORD nx1,nx2,nxx;
944 ttstop = t + t[1]; tt = t+FUNHEAD;
945 while ( tt < ttstop ) {
946 narg++;
947 if ( narg == 1 ) arg1 = tt;
948 else arg2 = tt;
949 NEXTARG(tt);
950 }
951 if ( narg != 2 ) goto defaultcase;
952 if ( *arg2 == -SNUMBER && arg2[1] <= 1 ) goto defaultcase;
953 else if ( *arg2 > 0 && ttstop[-1] < 0 ) goto defaultcase;
954 x1 = NumberMalloc("Norm-MakeRational");
955 if ( *arg1 == -SNUMBER ) {
956 if ( arg1[1] == 0 ) {
957 NumberFree(x1,"Norm-MakeRational");
958 goto NormZero;
959 }
960 if ( arg1[1] < 0 ) { x1[0] = -arg1[1]; nx1 = -1; }
961 else { x1[0] = arg1[1]; nx1 = 1; }
962 }
963 else if ( *arg1 > 0 ) {
964 WORD *tc;
965 nx1 = (ABS(arg2[-1])-1)/2;
966 tc = arg1+ARGHEAD+1+nx1;
967 if ( tc[0] != 1 ) {
968 NumberFree(x1,"Norm-MakeRational");
969 goto defaultcase;
970 }
971 for ( i = 1; i < nx1; i++ ) if ( tc[i] != 0 ) {
972 NumberFree(x1,"Norm-MakeRational");
973 goto defaultcase;
974 }
975 tc = arg1+ARGHEAD+1;
976 for ( i = 0; i < nx1; i++ ) x1[i] = tc[i];
977 if ( arg2[-1] < 0 ) nx1 = -nx1;
978 }
979 else {
980 NumberFree(x1,"Norm-MakeRational");
981 goto defaultcase;
982 }
983 x2 = NumberMalloc("Norm-MakeRational");
984 if ( *arg2 == -SNUMBER ) {
985 if ( arg2[1] <= 1 ) {
986 NumberFree(x2,"Norm-MakeRational");
987 NumberFree(x1,"Norm-MakeRational");
988 goto defaultcase;
989 }
990 else { x2[0] = arg2[1]; nx2 = 1; }
991 }
992 else if ( *arg2 > 0 ) {
993 WORD *tc;
994 nx2 = (ttstop[-1]-1)/2;
995 tc = arg2+ARGHEAD+1+nx2;
996 if ( tc[0] != 1 ) {
997 NumberFree(x2,"Norm-MakeRational");
998 NumberFree(x1,"Norm-MakeRational");
999 goto defaultcase;
1000 }
1001 for ( i = 1; i < nx2; i++ ) if ( tc[i] != 0 ) {
1002 NumberFree(x2,"Norm-MakeRational");
1003 NumberFree(x1,"Norm-MakeRational");
1004 goto defaultcase;
1005 }
1006 tc = arg2+ARGHEAD+1;
1007 for ( i = 0; i < nx2; i++ ) x2[i] = tc[i];
1008 }
1009 else {
1010 NumberFree(x2,"Norm-MakeRational");
1011 NumberFree(x1,"Norm-MakeRational");
1012 goto defaultcase;
1013 }
1014 if ( BigLong(x1,ABS(nx1),x2,nx2) >= 0 ) {
1015 UWORD *x3 = NumberMalloc("Norm-MakeRational");
1016 UWORD *x4 = NumberMalloc("Norm-MakeRational");
1017 WORD nx3, nx4;
1018 DivLong(x1,nx1,x2,nx2,x3,&nx3,x4,&nx4);
1019 for ( i = 0; i < ABS(nx4); i++ ) x1[i] = x4[i];
1020 nx1 = nx4;
1021 NumberFree(x4,"Norm-MakeRational");
1022 NumberFree(x3,"Norm-MakeRational");
1023 }
1024 xx = (UWORD *)(TermMalloc("Norm-MakeRational"));
1025 if ( MakeLongRational(BHEAD x1,nx1,x2,nx2,xx,&nxx) ) {
1026 static int warnflag = 1;
1027 if ( warnflag ) {
1028 MesPrint("%w Warning: fraction could not be reconstructed in MakeRational_");
1029 warnflag = 0;
1030 }
1031 ncoef = REDLENG(ncoef);
1032 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,x1,nx1) )
1033 goto FromNorm;
1034 }
1035 else {
1036 ncoef = REDLENG(ncoef);
1037 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,xx,nxx,
1038 (UWORD *)n_coef,&ncoef) ) goto FromNorm;
1039 }
1040 ncoef = INCLENG(ncoef);
1041 TermFree(xx,"Norm-MakeRational");
1042 NumberFree(x2,"Norm-MakeRational");
1043 NumberFree(x1,"Norm-MakeRational");
1044 }
1045 break;
1046 case TERMFUNCTION:
1047 if ( t[1] == FUNHEAD && AN.cTerm ) {
1048 ANsr = r; ANsm = m; ANsc = AN.cTerm;
1049 AN.cTerm = 0;
1050 t = ANsc + 1;
1051 m = ANsc + *ANsc;
1052 ncoef = REDLENG(ncoef);
1053 nnum = REDLENG(m[-1]);
1054 m -= ABS(m[-1]);
1055 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,nnum,
1056 (UWORD *)n_coef,&ncoef) ) goto FromNorm;
1057 ncoef = INCLENG(ncoef);
1058 r = t;
1059 }
1060 break;
1061 case FIRSTBRACKET:
1062 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1063 if ( GetFirstBracket(termout,t[FUNHEAD+1]) < 0 ) goto FromNorm;
1064 if ( *termout == 0 ) goto NormZero;
1065 if ( *termout > 4 ) {
1066 WORD *r1, *r2, *r3;
1067 while ( r < m ) *t++ = *r++;
1068 r1 = term + *term;
1069 r2 = termout + *termout; r2 -= ABS(r2[-1]);
1070 while ( r < r1 ) *r2++ = *r++;
1071 r3 = termout + 1;
1072 while ( r3 < r2 ) *t++ = *r3++;
1073 *term = t - term;
1074 if ( AT.WorkPointer > term && AT.WorkPointer < t )
1075 AT.WorkPointer = t;
1076 goto Restart;
1077 }
1078 }
1079 break;
1080 case FIRSTTERM:
1081 case CONTENTTERM:
1082 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1083 {
1084 EXPRESSIONS e = Expressions+t[FUNHEAD+1];
1085 POSITION oldondisk = AS.OldOnFile[t[FUNHEAD+1]];
1086 if ( e->replace == NEWLYDEFINEDEXPRESSION ) {
1087 AS.OldOnFile[t[FUNHEAD+1]] = e->onfile;
1088 }
1089 if ( *t == FIRSTTERM ) {
1090 if ( GetFirstTerm(termout,t[FUNHEAD+1],0) < 0 ) goto FromNorm;
1091 }
1092 else if ( *t == CONTENTTERM ) {
1093 if ( GetContent(termout,t[FUNHEAD+1]) < 0 ) goto FromNorm;
1094 }
1095 AS.OldOnFile[t[FUNHEAD+1]] = oldondisk;
1096 if ( *termout == 0 ) goto NormZero;
1097 }
1098PasteIn:;
1099 {
1100 WORD *r1, *r2, *r3, *r4, *r5, nr1, *rterm;
1101 r2 = termout + *termout; lnum = r2 - ABS(r2[-1]);
1102 nnum = REDLENG(r2[-1]);
1103
1104 r1 = term + *term; r3 = r1 - ABS(r1[-1]);
1105 nr1 = REDLENG(r1[-1]);
1106 if ( Mully(BHEAD (UWORD *)lnum,&nnum,(UWORD *)r3,nr1) ) goto FromNorm;
1107 nnum = INCLENG(nnum); nr1 = ABS(nnum); lnum[nr1-1] = nnum;
1108 rterm = TermMalloc("FirstTerm/ContentTerm");
1109 r4 = rterm+1; r5 = term+1; while ( r5 < t ) *r4++ = *r5++;
1110 r5 = termout+1; while ( r5 < lnum ) *r4++ = *r5++;
1111 r5 = r; while ( r5 < r3 ) *r4++ = *r5++;
1112 r5 = lnum; NCOPY(r4,r5,nr1);
1113 *rterm = r4-rterm;
1114 nr1 = *rterm; r1 = term; r2 = rterm; NCOPY(r1,r2,nr1);
1115 TermFree(rterm,"FirstTerm/ContentTerm");
1116 if ( AT.WorkPointer > term && AT.WorkPointer < r1 )
1117 AT.WorkPointer = r1;
1118 goto Restart;
1119 }
1120 }
1121 else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
1122 DOLLARS d = Dollars + t[FUNHEAD+1], newd = 0;
1123 int idol, ido;
1124#ifdef WITHPTHREADS
1125 int nummodopt, dtype = -1;
1126 if ( AS.MultiThreaded && ( AC.mparallelflag == PARALLELFLAG ) ) {
1127 for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
1128 if ( t[FUNHEAD+1] == ModOptdollars[nummodopt].number ) break;
1129 }
1130 if ( nummodopt < NumModOptdollars ) {
1131 dtype = ModOptdollars[nummodopt].type;
1132 if ( DollarLocalCopy(dtype) ) {
1133 d = ModOptdollars[nummodopt].dstruct+AT.identity;
1134 }
1135 }
1136 }
1137#endif
1138 if ( d->where && ( d->type == DOLTERMS || d->type == DOLNUMBER ) ) {
1139 newd = d;
1140 }
1141 else {
1142 if ( ( newd = DolToTerms(BHEAD t[FUNHEAD+1]) ) == 0 )
1143 goto NormZero;
1144 }
1145 if ( newd->where[0] == 0 ) {
1146 M_free(newd,"Copy of dollar variable");
1147 goto NormZero;
1148 }
1149 if ( *t == FIRSTTERM ) {
1150 idol = newd->where[0];
1151 for ( ido = 0; ido < idol; ido++ ) termout[ido] = newd->where[ido];
1152 }
1153 else if ( *t == CONTENTTERM ) {
1154 WORD *tterm;
1155 tterm = newd->where;
1156 idol = tterm[0];
1157 for ( ido = 0; ido < idol; ido++ ) termout[ido] = tterm[ido];
1158 tterm += *tterm;
1159 while ( *tterm ) {
1160 if ( ContentMerge(BHEAD termout,tterm) < 0 ) goto FromNorm;
1161 tterm += *tterm;
1162 }
1163 }
1164 if ( newd != d ) {
1165 if ( newd->factors ) M_free(newd->factors,"Dollar factors");
1166 M_free(newd,"Copy of dollar variable");
1167 newd = 0;
1168 }
1169 goto PasteIn;
1170 }
1171 break;
1172 case TERMSINEXPR:
1173 {
1174 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1175 x = TermsInExpression(t[FUNHEAD+1]);
1176multermnum: if ( x == 0 ) goto NormZero;
1177 if ( x < 0 ) {
1178 x = -x;
1179 if ( x > (LONG)WORDMASK ) { lnum[0] = x & WORDMASK;
1180 lnum[1] = x >> BITSINWORD; nnum = -2;
1181 }
1182 else { lnum[0] = x; nnum = -1; }
1183 }
1184 else if ( x > (LONG)WORDMASK ) {
1185 lnum[0] = x & WORDMASK;
1186 lnum[1] = x >> BITSINWORD;
1187 nnum = 2;
1188 }
1189 else { lnum[0] = x; nnum = 1; }
1190 ncoef = REDLENG(ncoef);
1191 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1192 goto FromNorm;
1193 ncoef = INCLENG(ncoef);
1194 }
1195 else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
1196 x = TermsInDollar(t[FUNHEAD+1]);
1197 goto multermnum;
1198 }
1199 else { pcom[ncom++] = t; }
1200 }
1201 break;
1202 case SIZEOFFUNCTION:
1203 {
1204 if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -EXPRESSION ) {
1205 x = SizeOfExpression(t[FUNHEAD+1]);
1206 goto multermnum;
1207 }
1208 else if ( ( t[1] == FUNHEAD+2 ) && t[FUNHEAD] == -DOLLAREXPRESSION ) {
1209 x = SizeOfDollar(t[FUNHEAD+1]);
1210 goto multermnum;
1211 }
1212 else { pcom[ncom++] = t; }
1213 }
1214 break;
1215 case MATCHFUNCTION:
1216 case PATTERNFUNCTION:
1217 break;
1218 case BINOMIAL:
1219/*
1220 Binomial function for internal use for the moment.
1221 The routine in reken.c should be more efficient.
1222*/
1223 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1224 && t[FUNHEAD+1] >= 0 && t[FUNHEAD+2] == -SNUMBER
1225 && t[FUNHEAD+3] >= 0 && t[FUNHEAD+1] >= t[FUNHEAD+3] ) {
1226 if ( t[FUNHEAD+1] > t[FUNHEAD+3] ) {
1227 if ( GetBinom((UWORD *)lnum,&nnum,
1228 t[FUNHEAD+1],t[FUNHEAD+3]) ) goto FromNorm;
1229 if ( nnum == 0 ) goto NormZero;
1230 goto MulIn;
1231 }
1232 }
1233 else pcom[ncom++] = t;
1234 break;
1235 case SIGNFUN:
1236/*
1237 Numerical function giving (-1)^arg
1238*/
1239 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1240 if ( ( t[FUNHEAD+1] & 1 ) != 0 ) ncoef = -ncoef;
1241 }
1242 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] )
1243 && ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) ) {
1244 UWORD *numer1,*denom1;
1245 WORD nsize = abs(t[t[1]-1]), nnsize, isize;
1246 nnsize = (nsize-1)/2;
1247 numer1 = (UWORD *)(t + FUNHEAD+ARGHEAD+1);
1248 denom1 = numer1 + nnsize;
1249 for ( isize = 1; isize < nnsize; isize++ ) {
1250 if ( denom1[isize] ) break;
1251 }
1252 if ( ( denom1[0] != 1 ) || isize < nnsize ) {
1253 pcom[ncom++] = t;
1254 }
1255 else {
1256 if ( ( numer1[0] & 1 ) != 0 ) ncoef = -ncoef;
1257 }
1258 }
1259 else {
1260 goto doflags;
1261/* pcom[ncom++] = t; */
1262 }
1263 break;
1264 case SIGFUNCTION:
1265/*
1266 Numerical function giving the sign of the numerical argument
1267 The sign of zero is 1.
1268 If there are roots of unity they are part of the sign.
1269*/
1270 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1271 if ( t[FUNHEAD+1] < 0 ) ncoef = -ncoef;
1272 }
1273 else if ( ( t[1] == FUNHEAD+2 ) && ( t[FUNHEAD] == -SYMBOL )
1274 && ( ( t[FUNHEAD+1] <= NumSymbols && t[FUNHEAD+1] > -MAXPOWER )
1275 && ( symbols[t[FUNHEAD+1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) ) {
1276 k = t[FUNHEAD+1];
1277 from = m;
1278 i = nsym;
1279 m = ppsym;
1280 if ( i > 0 ) do {
1281 m -= 2;
1282 if ( k == *m ) {
1283 m++;
1284 *m = *m + 1;
1285 *m %= symbols[k].maxpower;
1286 if ( ( symbols[k].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
1287 if ( ( ( symbols[k].maxpower & 1 ) == 0 ) &&
1288 ( *m >= symbols[k].maxpower/2 ) ) {
1289 *m -= symbols[k].maxpower/2; ncoef = -ncoef;
1290 }
1291 }
1292 if ( !*m ) {
1293 m--;
1294 while ( i < nsym )
1295 { *m = m[2]; m++; *m = m[2]; m++; i++; }
1296 ppsym -= 2;
1297 nsym--;
1298 }
1299 goto sigDoneSymbol;
1300 }
1301 } while ( k < *m && --i > 0 );
1302 m = ppsym;
1303 while ( i < nsym )
1304 { m--; m[2] = *m; m--; m[2] = *m; i++; }
1305 *m++ = k;
1306 *m = 1;
1307 ppsym += 2;
1308 nsym++;
1309sigDoneSymbol:;
1310 m = from;
1311 }
1312 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] ) ) {
1313 if ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) {
1314 if ( t[t[1]-1] < 0 ) ncoef = -ncoef;
1315 }
1316/*
1317 Now we should fish out the roots of unity
1318*/
1319 else if ( ( t[FUNHEAD+ARGHEAD]+FUNHEAD+ARGHEAD == t[1] )
1320 && ( t[FUNHEAD+ARGHEAD+1] == SYMBOL ) ) {
1321 WORD *ts = t + FUNHEAD+ARGHEAD+3;
1322 WORD its = ts[-1]-2;
1323 while ( its > 0 ) {
1324 if ( ( *ts != 0 ) && ( ( *ts > NumSymbols || *ts <= -MAXPOWER )
1325 || ( symbols[*ts].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY ) ) {
1326 goto signogood;
1327 }
1328 ts += 2; its -= 2;
1329 }
1330/*
1331 Now we have only roots of unity which should be
1332 registered in the list of symbols.
1333*/
1334 if ( t[t[1]-1] < 0 ) ncoef = -ncoef;
1335 ts = t + FUNHEAD+ARGHEAD+3;
1336 its = ts[-1]-2;
1337 from = m;
1338 while ( its > 0 ) {
1339 i = nsym;
1340 m = ppsym;
1341 if ( i > 0 ) do {
1342 m -= 2;
1343 if ( *ts == *m ) {
1344 ts++; m++;
1345 *m += *ts;
1346 if ( ( ts[-1] <= NumSymbols && ts[-1] > -MAXPOWER ) &&
1347 ( symbols[ts[-1]].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
1348 *m %= symbols[ts[-1]].maxpower;
1349 if ( *m < 0 ) *m += symbols[ts[-1]].maxpower;
1350 if ( ( symbols[ts[-1]].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
1351 if ( ( ( symbols[ts[-1]].maxpower & 1 ) == 0 ) &&
1352 ( *m >= symbols[ts[-1]].maxpower/2 ) ) {
1353 *m -= symbols[ts[-1]].maxpower/2; ncoef = -ncoef;
1354 }
1355 }
1356 }
1357 if ( !*m ) {
1358 m--;
1359 while ( i < nsym )
1360 { *m = m[2]; m++; *m = m[2]; m++; i++; }
1361 ppsym -= 2;
1362 nsym--;
1363 }
1364 ts++; its -= 2;
1365 goto sigNextSymbol;
1366 }
1367 } while ( *ts < *m && --i > 0 );
1368 m = ppsym;
1369 while ( i < nsym )
1370 { m--; m[2] = *m; m--; m[2] = *m; i++; }
1371 *m++ = *ts++;
1372 *m = *ts++;
1373 ppsym += 2;
1374 nsym++; its -= 2;
1375sigNextSymbol:;
1376 }
1377 m = from;
1378 }
1379 else {
1380signogood: pcom[ncom++] = t;
1381 }
1382 }
1383 else pcom[ncom++] = t;
1384 break;
1385 case ABSFUNCTION:
1386/*
1387 Numerical function giving the absolute value of the
1388 numerical argument. Or roots of unity.
1389*/
1390 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1391 k = t[FUNHEAD+1];
1392 if ( k < 0 ) k = -k;
1393 if ( k == 0 ) goto NormZero;
1394 *((UWORD *)lnum) = k; nnum = 1;
1395 goto MulIn;
1396
1397 }
1398 else if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SYMBOL ) {
1399 k = t[FUNHEAD+1];
1400 if ( ( k > NumSymbols || k <= -MAXPOWER )
1401 || ( symbols[k].complex & VARTYPEROOTOFUNITY ) != VARTYPEROOTOFUNITY )
1402 goto absnogood;
1403 }
1404 else if ( ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+t[FUNHEAD] )
1405 && ( t[1] == FUNHEAD+ARGHEAD+t[FUNHEAD+ARGHEAD] ) ) {
1406 if ( t[FUNHEAD] == ARGHEAD+1+abs(t[t[1]-1]) ) {
1407 WORD *ts;
1408absnosymbols: ts = t + t[1] -1;
1409 ncoef = REDLENG(ncoef);
1410 nnum = REDLENG(*ts);
1411 if ( nnum < 0 ) nnum = -nnum;
1412 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,
1413 (UWORD *)(ts-ABS(*ts)+1),nnum,
1414 (UWORD *)n_coef,&ncoef) ) goto FromNorm;
1415 ncoef = INCLENG(ncoef);
1416 }
1417/*
1418 Now get rid of the roots of unity. This includes i_
1419*/
1420 else if ( t[FUNHEAD+ARGHEAD+1] == SYMBOL ) {
1421 WORD *ts = t+FUNHEAD+ARGHEAD+1;
1422 WORD its = ts[1] - 2;
1423 ts += 2;
1424 while ( its > 0 ) {
1425 if ( *ts == 0 ) { }
1426 else if ( ( *ts > NumSymbols || *ts <= -MAXPOWER )
1427 || ( symbols[*ts].complex & VARTYPEROOTOFUNITY )
1428 != VARTYPEROOTOFUNITY ) goto absnogood;
1429 its -= 2; ts += 2;
1430 }
1431 goto absnosymbols;
1432 }
1433 else {
1434absnogood: pcom[ncom++] = t;
1435 }
1436 }
1437 else pcom[ncom++] = t;
1438 break;
1439 case MODFUNCTION:
1440 case MOD2FUNCTION:
1441/*
1442 Mod function. Does work if two arguments and the
1443 second argument is a positive short number
1444*/
1445 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1446 && t[FUNHEAD+2] == -SNUMBER && t[FUNHEAD+3] > 1 ) {
1447 WORD tmod;
1448 tmod = (t[FUNHEAD+1]%t[FUNHEAD+3]);
1449 if ( tmod < 0 ) tmod += t[FUNHEAD+3];
1450 if ( *t == MOD2FUNCTION && tmod > t[FUNHEAD+3]/2 )
1451 tmod -= t[FUNHEAD+3];
1452 if ( tmod < 0 ) {
1453 *((UWORD *)lnum) = -tmod;
1454 nnum = -1;
1455 }
1456 else if ( tmod > 0 ) {
1457 *((UWORD *)lnum) = tmod;
1458 nnum = 1;
1459 }
1460 else goto NormZero;
1461 goto MulIn;
1462 }
1463 else if ( t[1] > t[FUNHEAD+2] && t[FUNHEAD] > 0
1464 && t[FUNHEAD+t[FUNHEAD]] == -SNUMBER
1465 && t[FUNHEAD+t[FUNHEAD]+1] > 1
1466 && t[1] == FUNHEAD+2+t[FUNHEAD] ) {
1467 WORD *ttt = t+FUNHEAD, iii;
1468 iii = ttt[*ttt-1];
1469 if ( *ttt == ttt[ARGHEAD]+ARGHEAD &&
1470 ttt[ARGHEAD] == ABS(iii)+1 ) {
1471 WORD ncmod = 1;
1472 WORD cmod = ttt[*ttt+1];
1473 iii = REDLENG(iii);
1474 if ( *t == MODFUNCTION ) {
1475 if ( TakeModulus((UWORD *)(ttt+ARGHEAD+1)
1476 ,&iii,(UWORD *)(&cmod),ncmod,UNPACK|NOINVERSES) )
1477 goto FromNorm;
1478 }
1479 else {
1480 if ( TakeModulus((UWORD *)(ttt+ARGHEAD+1)
1481 ,&iii,(UWORD *)(&cmod),ncmod,UNPACK|POSNEG|NOINVERSES) )
1482 goto FromNorm;
1483 }
1484 if ( *t == MOD2FUNCTION && ttt[ARGHEAD+1] > cmod/2 && iii > 0 ) {
1485 ttt[ARGHEAD+1] -= cmod;
1486 }
1487 if ( ttt[ARGHEAD+1] < 0 ) {
1488 *((UWORD *)lnum) = -ttt[ARGHEAD+1];
1489 nnum = -1;
1490 }
1491 else if ( ttt[ARGHEAD+1] > 0 ) {
1492 *((UWORD *)lnum) = ttt[ARGHEAD+1];
1493 nnum = 1;
1494 }
1495 else goto NormZero;
1496 goto MulIn;
1497 }
1498 }
1499 else if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER ) {
1500 *((UWORD *)lnum) = t[FUNHEAD+1];
1501 if ( *lnum == 0 ) goto NormZero;
1502 nnum = 1;
1503 goto MulIn;
1504 }
1505 else if ( ( ( t[FUNHEAD] < 0 ) && ( t[FUNHEAD] == -SNUMBER )
1506 && ( t[1] >= ( FUNHEAD+6+ARGHEAD ) )
1507 && ( t[FUNHEAD+2] >= 4+ARGHEAD )
1508 && ( t[t[1]-1] == t[FUNHEAD+2+ARGHEAD]-1 ) ) ||
1509 ( ( t[FUNHEAD] > 0 )
1510 && ( t[FUNHEAD]-ARGHEAD-1 == ABS(t[FUNHEAD+t[FUNHEAD]-1]) )
1511 && ( t[FUNHEAD+t[FUNHEAD]]-ARGHEAD-1 == t[t[1]-1] ) ) ) {
1512/*
1513 Check that the last (long) number is integer
1514*/
1515 WORD *ttt = t + t[1], iii, iii1;
1516 UWORD coefbuf[2], *coef2, ncoef2;
1517 iii = (ttt[-1]-1)/2;
1518 ttt -= iii;
1519 if ( ttt[-1] != 1 ) {
1520exitfromhere:
1521 pcom[ncom++] = t;
1522 break;
1523 }
1524 iii--;
1525 for ( iii1 = 0; iii1 < iii; iii1++ ) {
1526 if ( ttt[iii1] != 0 ) goto exitfromhere;
1527 }
1528/*
1529 Now we have a hit!
1530 The first argument will be put in lnum.
1531 It will be a rational.
1532 The second argument will be a long integer in coef2.
1533*/
1534 ttt = t + FUNHEAD;
1535 if ( *ttt < 0 ) {
1536 if ( ttt[1] < 0 ) {
1537 nnum = -1; lnum[0] = -ttt[1]; lnum[1] = 1;
1538 }
1539 else {
1540 nnum = 1; lnum[0] = ttt[1]; lnum[1] = 1;
1541 }
1542 }
1543 else {
1544 nnum = ABS(ttt[ttt[0]-1] - 1);
1545 for ( iii = 0; iii < nnum; iii++ ) {
1546 lnum[iii] = ttt[ARGHEAD+1+iii];
1547 }
1548 nnum = nnum/2;
1549 if ( ttt[ttt[0]-1] < 0 ) nnum = -nnum;
1550 }
1551 NEXTARG(ttt);
1552 if ( *ttt < 0 ) {
1553 coef2 = coefbuf;
1554 ncoef2 = 3; *coef2 = (UWORD)(ttt[1]);
1555 coef2[1] = 1;
1556 }
1557 else {
1558 coef2 = (UWORD *)(ttt+ARGHEAD+1);
1559 ncoef2 = (ttt[ttt[0]-1]-1)/2;
1560 }
1561 if ( TakeModulus((UWORD *)lnum,&nnum,(UWORD *)coef2,ncoef2,
1562 UNPACK|NOINVERSES|FROMFUNCTION) ) {
1563 goto FromNorm;
1564 }
1565 if ( *t == MOD2FUNCTION && nnum > 0 ) {
1566 UWORD *coef3 = NumberMalloc("Mod2Function"), two = 2;
1567 WORD ncoef3;
1568 if ( MulLong((UWORD *)lnum,nnum,&two,1,coef3,&ncoef3) )
1569 goto FromNorm;
1570 if ( BigLong(coef3,ncoef3,(UWORD *)coef2,ncoef2) > 0 ) {
1571 nnum = -nnum;
1572 AddLong((UWORD *)lnum,nnum,(UWORD *)coef2,ncoef2
1573 ,(UWORD *)lnum,&nnum);
1574 nnum = -nnum;
1575 }
1576 NumberFree(coef3,"Mod2Function");
1577 }
1578/*
1579 Do we have to pack? No, because the answer is not a fraction
1580*/
1581 ncoef = REDLENG(ncoef);
1582 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
1583 goto FromNorm;
1584 ncoef = INCLENG(ncoef);
1585 }
1586 else pcom[ncom++] = t;
1587 break;
1588 case EXTEUCLIDEAN:
1589 {
1590 WORD argcount = 0, *tc, *ts, xc, xs, *tcc;
1591 UWORD *Num1, *Num2, *Num3, *Num4;
1592 WORD size1, size2, size3, size4, space;
1593 tc = t+FUNHEAD; ts = t + t[1];
1594 while ( argcount < 3 && tc < ts ) { NEXTARG(tc); argcount++; }
1595 if ( argcount != 2 ) goto defaultcase;
1596 if ( t[FUNHEAD] == -SNUMBER ) {
1597 if ( t[FUNHEAD+1] <= 1 ) goto defaultcase;
1598 if ( t[FUNHEAD+2] == -SNUMBER ) {
1599 if ( t[FUNHEAD+3] <= 1 ) goto defaultcase;
1600 Num2 = NumberMalloc("modinverses");
1601 *Num2 = t[FUNHEAD+3]; size2 = 1;
1602 }
1603 else {
1604 if ( ts[-1] < 0 ) goto defaultcase;
1605 if ( ts[-1] != t[FUNHEAD+2]-ARGHEAD-1 ) goto defaultcase;
1606 xs = (ts[-1]-1)/2;
1607 tcc = ts-xs-1;
1608 if ( *tcc != 1 ) goto defaultcase;
1609 for ( i = 1; i < xs; i++ ) {
1610 if ( tcc[i] != 0 ) goto defaultcase;
1611 }
1612 Num2 = NumberMalloc("modinverses");
1613 size2 = xs;
1614 for ( i = 0; i < xs; i++ ) Num2[i] = t[FUNHEAD+ARGHEAD+3+i];
1615 }
1616 Num1 = NumberMalloc("modinverses");
1617 *Num1 = t[FUNHEAD+1]; size1 = 1;
1618 }
1619 else {
1620 tc = t + FUNHEAD + t[FUNHEAD];
1621 if ( tc[-1] < 0 ) goto defaultcase;
1622 if ( tc[-1] != t[FUNHEAD]-ARGHEAD-1 ) goto defaultcase;
1623 xc = (tc[-1]-1)/2;
1624 tcc = tc-xc-1;
1625 if ( *tcc != 1 ) goto defaultcase;
1626 for ( i = 1; i < xc; i++ ) {
1627 if ( tcc[i] != 0 ) goto defaultcase;
1628 }
1629 if ( *tc == -SNUMBER ) {
1630 if ( tc[1] <= 1 ) goto defaultcase;
1631 Num2 = NumberMalloc("modinverses");
1632 *Num2 = tc[1]; size2 = 1;
1633 }
1634 else {
1635 if ( ts[-1] < 0 ) goto defaultcase;
1636 if ( ts[-1] != t[FUNHEAD+2]-ARGHEAD-1 ) goto defaultcase;
1637 xs = (ts[-1]-1)/2;
1638 tcc = ts-xs-1;
1639 if ( *tcc != 1 ) goto defaultcase;
1640 for ( i = 1; i < xs; i++ ) {
1641 if ( tcc[i] != 0 ) goto defaultcase;
1642 }
1643 Num2 = NumberMalloc("modinverses");
1644 size2 = xs;
1645 for ( i = 0; i < xs; i++ ) Num2[i] = tc[ARGHEAD+1+i];
1646 }
1647 Num1 = NumberMalloc("modinverses");
1648 size1 = xc;
1649 for ( i = 0; i < xc; i++ ) Num1[i] = t[FUNHEAD+ARGHEAD+1+i];
1650 }
1651 Num3 = NumberMalloc("modinverses");
1652 Num4 = NumberMalloc("modinverses");
1653 GetLongModInverses(BHEAD Num1,size1,Num2,size2
1654 ,Num3,&size3,Num4,&size4);
1655/*
1656 Now we have to compose the answer. This needs more space
1657 and hence we have to put this inside the term.
1658 Compute first how much extra space we need.
1659 Then move the trailing part of the term upwards.
1660 Do not forget relevant pointers!!! (r, m, termout, AT.WorkPointer)
1661*/
1662 space = 0;
1663 if ( ( size3 == 1 || size3 == -1 ) && (*Num3&TOPBITONLY) == 0 ) space += 2;
1664 else space += ARGHEAD + 2*ABS(size3) + 2;
1665 if ( ( size4 == 1 || size4 == -1 ) && (*Num4&TOPBITONLY) == 0 ) space += 2;
1666 else space += ARGHEAD + 2*ABS(size4) + 2;
1667 tt = term + *term; u = tt + space;
1668 while ( tt >= ts ) *--u = *--tt;
1669 m += space; r += space;
1670 *term += space;
1671 t[1] += space;
1672 if ( ( size3 == 1 || size3 == -1 ) && (*Num3&TOPBITONLY) == 0 ) {
1673 *ts++ = -SNUMBER; *ts = (WORD)(*Num3);
1674 if ( size3 < 0 ) *ts = -*ts;
1675 ts++;
1676 }
1677 else {
1678 *ts++ = 2*ABS(size3)+ARGHEAD+2;
1679 *ts++ = 0; FILLARG(ts)
1680 *ts++ = 2*ABS(size3)+1;
1681 for ( i = 0; i < ABS(size3); i++ ) *ts++ = Num3[i];
1682 *ts++ = 1;
1683 for ( i = 1; i < ABS(size3); i++ ) *ts++ = 0;
1684 if ( size3 < 0 ) *ts++ = 2*size3-1;
1685 else *ts++ = 2*size3+1;
1686 }
1687 if ( ( size4 == 1 || size4 == -1 ) && (*Num4&TOPBITONLY) == 0 ) {
1688 *ts++ = -SNUMBER; *ts = *Num4;
1689 if ( size4 < 0 ) *ts = -*ts;
1690 ts++;
1691 }
1692 else {
1693 *ts++ = 2*ABS(size4)+ARGHEAD+2;
1694 *ts++ = 0; FILLARG(ts)
1695 *ts++ = 2*ABS(size4)+2;
1696 for ( i = 0; i < ABS(size4); i++ ) *ts++ = Num4[i];
1697 *ts++ = 1;
1698 for ( i = 1; i < ABS(size4); i++ ) *ts++ = 0;
1699 if ( size4 < 0 ) *ts++ = 2*size4-1;
1700 else *ts++ = 2*size4+1;
1701 }
1702 NumberFree(Num4,"modinverses");
1703 NumberFree(Num3,"modinverses");
1704 NumberFree(Num1,"modinverses");
1705 NumberFree(Num2,"modinverses");
1706 t[2] = 0; /* mark function as clean. */
1707 goto Restart;
1708 }
1709 break;
1710 case GCDFUNCTION:
1711#ifdef EVALUATEGCD
1712#ifdef NEWGCDFUNCTION
1713 {
1714/*
1715 Has two integer arguments
1716 Four cases: S,S, S,L, L,S, L,L
1717*/
1718 WORD *num1, *num2, size1, size2, stor1, stor2, *ttt, ti;
1719 if ( t[1] == FUNHEAD+4 && t[FUNHEAD] == -SNUMBER
1720 && t[FUNHEAD+2] == -SNUMBER && t[FUNHEAD+1] != 0
1721 && t[FUNHEAD+3] != 0 ) { /* Short,Short */
1722 stor1 = t[FUNHEAD+1];
1723 stor2 = t[FUNHEAD+3];
1724 if ( stor1 < 0 ) stor1 = -stor1;
1725 if ( stor2 < 0 ) stor2 = -stor2;
1726 num1 = &stor1; num2 = &stor2;
1727 size1 = size2 = 1;
1728 goto gcdcalc;
1729 }
1730 else if ( t[1] > FUNHEAD+4 ) {
1731 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] != 0
1732 && t[FUNHEAD+2] == t[1]-FUNHEAD-2 &&
1733 ABS(t[t[1]-1]) == t[FUNHEAD+2]-1-ARGHEAD ) { /* Short,Long */
1734 num2 = t + t[1];
1735 size2 = ABS(num2[-1]);
1736 ttt = num2-1;
1737 num2 -= size2;
1738 size2 = (size2-1)/2;
1739 ti = size2;
1740 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1741 if ( ti == 1 && ttt[-1] == 1 ) {
1742 stor1 = t[FUNHEAD+1];
1743 if ( stor1 < 0 ) stor1 = -stor1;
1744 num1 = &stor1;
1745 size1 = 1;
1746 goto gcdcalc;
1747 }
1748 else pcom[ncom++] = t;
1749 }
1750 else if ( t[FUNHEAD] > 0 &&
1751 t[FUNHEAD]-1-ARGHEAD == ABS(t[t[FUNHEAD]+FUNHEAD-1]) ) {
1752 num1 = t + FUNHEAD + t[FUNHEAD];
1753 size1 = ABS(num1[-1]);
1754 ttt = num1-1;
1755 num1 -= size1;
1756 size1 = (size1-1)/2;
1757 ti = size1;
1758 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1759 if ( ti == 1 && ttt[-1] == 1 ) {
1760 if ( t[1]-FUNHEAD == t[FUNHEAD]+2 && t[t[1]-2] == -SNUMBER
1761 && t[t[1]-1] != 0 ) { /* Long,Short */
1762 stor2 = t[t[1]-1];
1763 if ( stor2 < 0 ) stor2 = -stor2;
1764 num2 = &stor2;
1765 size2 = 1;
1766 goto gcdcalc;
1767 }
1768 else if ( t[1]-FUNHEAD == t[FUNHEAD]+t[FUNHEAD+t[FUNHEAD]]
1769 && ABS(t[t[1]-1]) == t[FUNHEAD+t[FUNHEAD]] - ARGHEAD-1 ) {
1770 num2 = t + t[1];
1771 size2 = ABS(num2[-1]);
1772 ttt = num2-1;
1773 num2 -= size2;
1774 size2 = (size2-1)/2;
1775 ti = size2;
1776 while ( ti > 1 && ttt[-1] == 0 ) { ttt--; ti--; }
1777 if ( ti == 1 && ttt[-1] == 1 ) {
1778gcdcalc: if ( GcdLong(BHEAD (UWORD *)num1,size1,(UWORD *)num2,size2
1779 ,(UWORD *)lnum,&nnum) ) goto FromNorm;
1780 goto MulIn;
1781 }
1782 else pcom[ncom++] = t;
1783 }
1784 else pcom[ncom++] = t;
1785 }
1786 else pcom[ncom++] = t;
1787 }
1788 else pcom[ncom++] = t;
1789 }
1790 else pcom[ncom++] = t;
1791 }
1792#else
1793 {
1794 WORD *gcd = AT.WorkPointer;
1795 if ( ( gcd = EvaluateGcd(BHEAD t) ) == 0 ) goto FromNorm;
1796 if ( *gcd == 4 && gcd[1] == 1 && gcd[2] == 1 && gcd[4] == 0 ) {
1797 AT.WorkPointer = gcd;
1798 }
1799 else if ( gcd[*gcd] == 0 ) {
1800 WORD *t1, iii, change, *num, *den, numsize, densize;
1801 if ( gcd[*gcd-1] < *gcd-1 ) {
1802 t1 = gcd+1;
1803 for ( iii = 2; iii < t1[1]; iii += 2 ) {
1804 change = ExtraSymbol(t1[iii],t1[iii+1],nsym,ppsym,&ncoef);
1805 nsym += change;
1806 ppsym += change * 2;
1807 }
1808 }
1809 t1 = gcd + *gcd;
1810 iii = t1[-1]; num = t1-iii; numsize = (iii-1)/2;
1811 den = num + numsize; densize = numsize;
1812 while ( numsize > 1 && num[numsize-1] == 0 ) numsize--;
1813 while ( densize > 1 && den[densize-1] == 0 ) densize--;
1814 if ( numsize > 1 || num[0] != 1 ) {
1815 ncoef = REDLENG(ncoef);
1816 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)num,numsize) ) goto FromNorm;
1817 ncoef = INCLENG(ncoef);
1818 }
1819 if ( densize > 1 || den[0] != 1 ) {
1820 ncoef = REDLENG(ncoef);
1821 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)den,densize) ) goto FromNorm;
1822 ncoef = INCLENG(ncoef);
1823 }
1824 AT.WorkPointer = gcd;
1825 }
1826 else { /* a whole expression */
1827/*
1828 Action: Put the expression in a compiler buffer.
1829 Insert a SUBEXPRESSION subterm
1830 Set the return value of the routine such that in
1831 Generator the term gets sent again to TestSub.
1832
1833 1: put in C (ebufnum)
1834 2: after that the WorkSpace is free again.
1835 3: insert the SUBEXPRESSION
1836 4: copy the top part of the term down
1837*/
1838 LONG size = AT.WorkPointer - gcd;
1839
1840 ss = AddRHS(AT.ebufnum,1);
1841 while ( (ss + size + 10) > C->Top ) ss = DoubleCbuffer(AT.ebufnum,ss,13);
1842 tt = gcd;
1843 NCOPY(ss,tt,size);
1844 C->rhs[C->numrhs+1] = ss;
1845 C->Pointer = ss;
1846
1847 t[0] = SUBEXPRESSION;
1848 t[1] = SUBEXPSIZE;
1849 t[2] = C->numrhs;
1850 t[3] = 1;
1851 t[4] = AT.ebufnum;
1852 t += 5;
1853 tt = term + *term;
1854 while ( r < tt ) *t++ = *r++;
1855 *term = t - term;
1856
1857 regval = 1;
1858 goto RegEnd;
1859 }
1860 }
1861#endif
1862#else
1863 MesPrint(" Unexpected call to EvaluateGCD");
1864 Terminate(-1);
1865#endif
1866 break;
1867 case MINFUNCTION:
1868 case MAXFUNCTION:
1869 if ( t[1] == FUNHEAD ) break;
1870 {
1871 WORD *ttt = t + FUNHEAD;
1872 WORD *tttstop = t + t[1];
1873 WORD tterm[4], iii;
1874 while ( ttt < tttstop ) {
1875 if ( *ttt > 0 ) {
1876 if ( ttt[ARGHEAD]-1 > ABS(ttt[*ttt-1]) ) goto nospec;
1877 ttt += *ttt;
1878 }
1879 else {
1880 if ( *ttt != -SNUMBER ) goto nospec;
1881 ttt += 2;
1882 }
1883 }
1884/*
1885 Function has only numerical arguments
1886 Pick up the first argument.
1887*/
1888 ttt = t + FUNHEAD;
1889 if ( *ttt > 0 ) {
1890loadnew1:
1891 for ( iii = 0; iii < ttt[ARGHEAD]; iii++ )
1892 n_llnum[iii] = ttt[ARGHEAD+iii];
1893 ttt += *ttt;
1894 }
1895 else {
1896loadnew2:
1897 if ( ttt[1] == 0 ) {
1898 n_llnum[0] = n_llnum[1] = n_llnum[2] = n_llnum[3] = 0;
1899 }
1900 else {
1901 n_llnum[0] = 4;
1902 if ( ttt[1] > 0 ) { n_llnum[1] = ttt[1]; n_llnum[3] = 3; }
1903 else { n_llnum[1] = -ttt[1]; n_llnum[3] = -3; }
1904 n_llnum[2] = 1;
1905 }
1906 ttt += 2;
1907 }
1908/*
1909 Now loop over the other arguments
1910*/
1911 while ( ttt < tttstop ) {
1912 if ( *ttt > 0 ) {
1913 if ( n_llnum[0] == 0 ) {
1914 if ( ( *t == MINFUNCTION && ttt[*ttt-1] < 0 )
1915 || ( *t == MAXFUNCTION && ttt[*ttt-1] > 0 ) )
1916 goto loadnew1;
1917 }
1918 else {
1919 ttt += ARGHEAD;
1920 iii = CompCoef(n_llnum,ttt);
1921 if ( ( iii > 0 && *t == MINFUNCTION )
1922 || ( iii < 0 && *t == MAXFUNCTION ) ) {
1923 for ( iii = 0; iii < ttt[0]; iii++ )
1924 n_llnum[iii] = ttt[iii];
1925 }
1926 }
1927 ttt += *ttt;
1928 }
1929 else {
1930 if ( n_llnum[0] == 0 ) {
1931 if ( ( *t == MINFUNCTION && ttt[1] < 0 )
1932 || ( *t == MAXFUNCTION && ttt[1] > 0 ) )
1933 goto loadnew2;
1934 }
1935 else if ( ttt[1] == 0 ) {
1936 if ( ( *t == MINFUNCTION && n_llnum[*n_llnum-1] > 0 )
1937 || ( *t == MAXFUNCTION && n_llnum[*n_llnum-1] < 0 ) ) {
1938 n_llnum[0] = 0;
1939 }
1940 }
1941 else {
1942 tterm[0] = 4; tterm[2] = 1;
1943 if ( ttt[1] < 0 ) { tterm[1] = -ttt[1]; tterm[3] = -3; }
1944 else { tterm[1] = ttt[1]; tterm[3] = 3; }
1945 iii = CompCoef(n_llnum,tterm);
1946 if ( ( iii > 0 && *t == MINFUNCTION )
1947 || ( iii < 0 && *t == MAXFUNCTION ) ) {
1948 for ( iii = 0; iii < 4; iii++ )
1949 n_llnum[iii] = tterm[iii];
1950 }
1951 }
1952 ttt += 2;
1953 }
1954 }
1955 if ( n_llnum[0] == 0 ) goto NormZero;
1956 ncoef = REDLENG(ncoef);
1957 nnum = REDLENG(n_llnum[*n_llnum-1]);
1958 if ( MulRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)lnum,nnum,
1959 (UWORD *)n_coef,&ncoef) ) goto FromNorm;
1960 ncoef = INCLENG(ncoef);
1961 }
1962 break;
1963 case INVERSEFACTORIAL:
1964 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] >= 0 ) {
1965 if ( Factorial(BHEAD t[FUNHEAD+1],(UWORD *)lnum,&nnum) )
1966 goto FromNorm;
1967 ncoef = REDLENG(ncoef);
1968 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) ) goto FromNorm;
1969 ncoef = INCLENG(ncoef);
1970 }
1971 else {
1972nospec: pcom[ncom++] = t;
1973 }
1974 break;
1975 case MAXPOWEROF:
1976 if ( ( t[FUNHEAD] == -SYMBOL )
1977 && ( t[FUNHEAD+1] > 0 ) && ( t[1] == FUNHEAD+2 ) ) {
1978 *((UWORD *)lnum) = symbols[t[FUNHEAD+1]].maxpower;
1979 nnum = 1;
1980 goto MulIn;
1981 }
1982 else { pcom[ncom++] = t; }
1983 break;
1984 case MINPOWEROF:
1985 if ( ( t[FUNHEAD] == -SYMBOL )
1986 && ( t[FUNHEAD] > 0 ) && ( t[1] == FUNHEAD+2 ) ) {
1987 *((UWORD *)lnum) = symbols[t[FUNHEAD+1]].minpower;
1988 nnum = 1;
1989 goto MulIn;
1990 }
1991 else { pcom[ncom++] = t; }
1992 break;
1993 case PRIMENUMBER :
1994 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER
1995 && t[FUNHEAD+1] > 0 ) {
1996 UWORD xp = (UWORD)(NextPrime(BHEAD t[FUNHEAD+1]));
1997 ncoef = REDLENG(ncoef);
1998 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,&xp,1) ) goto FromNorm;
1999 ncoef = INCLENG(ncoef);
2000 }
2001 else goto defaultcase;
2002 break;
2003 case MOEBIUS:
2004/*
2005 Numerical function giving -1,0,1 or no evaluation
2006*/
2007 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER
2008 && t[FUNHEAD+1] > 0 ) {
2009 WORD val = Moebius(BHEAD t[FUNHEAD+1]);
2010 if ( val == 0 ) goto NormZero;
2011 if ( val < 0 ) ncoef = -ncoef;
2012 }
2013 else {
2014 pcom[ncom++] = t;
2015 }
2016 break;
2017 case LNUMBER :
2018 ncoef = REDLENG(ncoef);
2019 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(t+3),t[2]) ) goto FromNorm;
2020 ncoef = INCLENG(ncoef);
2021 break;
2022 case SNUMBER :
2023 if ( t[2] < 0 ) {
2024 t[2] = -t[2];
2025 if ( t[3] & 1 ) ncoef = -ncoef;
2026 }
2027 else if ( t[2] == 0 ) {
2028 if ( t[3] < 0 ) goto NormInf;
2029 goto NormZero;
2030 }
2031 lnum[0] = t[2];
2032 nnum = 1;
2033 if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) ) goto FromNorm;
2034 ncoef = REDLENG(ncoef);
2035 if ( t[3] < 0 ) {
2036 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2037 goto FromNorm;
2038 }
2039 else if ( t[3] > 0 ) {
2040 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2041 goto FromNorm;
2042 }
2043 ncoef = INCLENG(ncoef);
2044 break;
2045 case GAMMA :
2046 case GAMMAI :
2047 case GAMMAFIVE :
2048 case GAMMASIX :
2049 case GAMMASEVEN :
2050 if ( t[1] == FUNHEAD ) {
2051 MLOCK(ErrorMessageLock);
2052 MesPrint("Gamma matrix without spin line encountered.");
2053 MUNLOCK(ErrorMessageLock);
2054 goto NormMin;
2055 }
2056 pnco[nnco++] = t;
2057 t += FUNHEAD+1;
2058 goto ScanCont;
2059 case LEVICIVITA :
2060 peps[neps++] = t;
2061 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG ) {
2062 t[2] &= ~DIRTYFLAG;
2063 t[2] |= DIRTYSYMFLAG;
2064 }
2065 t += FUNHEAD;
2066ScanCont: while ( t < r ) {
2067 if ( *t >= AM.OffsetIndex &&
2068 ( *t >= AM.DumInd || ( *t < AM.WilInd &&
2069 indices[*t-AM.OffsetIndex].dimension ) ) )
2070 pcon[ncon++] = t;
2071 t++;
2072 }
2073 break;
2074 case EXPONENT :
2075 {
2076 WORD *rr;
2077 k = 1;
2078 rr = t + FUNHEAD;
2079 if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) )
2080 k = 0;
2081 if ( *rr == -SNUMBER && rr[1] == 1 ) break;
2082 if ( *rr <= -FUNCTION ) k = *rr;
2083 NEXTARG(rr)
2084 if ( *rr == ARGHEAD || ( *rr == -SNUMBER && rr[1] == 0 ) ) {
2085 if ( k == 0 ) goto NormZZ;
2086 break;
2087 }
2088 if ( *rr == -SNUMBER && rr[1] > 0 && rr[1] < MAXPOWER && k < 0 ) {
2089 k = -k;
2090 if ( functions[k-FUNCTION].commute ) {
2091 for ( i = 0; i < rr[1]; i++ ) pnco[nnco++] = rr-1;
2092 }
2093 else {
2094 for ( i = 0; i < rr[1]; i++ ) pcom[ncom++] = rr-1;
2095 }
2096 break;
2097 }
2098 if ( k == 0 ) goto NormZero;
2099 if ( t[FUNHEAD] == -SYMBOL && *rr == -SNUMBER && t[1] == FUNHEAD+4 ) {
2100 if ( rr[1] < MAXPOWER ) {
2101 t[FUNHEAD+2] = t[FUNHEAD+1]; t += FUNHEAD+2;
2102 from = m;
2103 goto NextSymbol;
2104 }
2105 }
2106/*
2107 if ( ( t[FUNHEAD] > 0 && t[FUNHEAD+1] != 0 )
2108 || ( *rr > 0 && rr[1] != 0 ) ) {}
2109 else
2110*/
2111 t[2] &= ~DIRTYSYMFLAG;
2112
2113 pnco[nnco++] = t;
2114 }
2115 break;
2116 case DENOMINATOR :
2117 t[2] &= ~DIRTYSYMFLAG;
2118 pden[nden++] = t;
2119 pnco[nnco++] = t;
2120 break;
2121 case INDEX :
2122 t += 2;
2123 do {
2124 if ( *t == 0 || *t == AM.vectorzero ) goto NormZero;
2125 if ( *t > 0 && *t < AM.OffsetIndex ) {
2126 lnum[0] = *t++;
2127 nnum = 1;
2128 ncoef = REDLENG(ncoef);
2129 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)lnum,nnum) )
2130 goto FromNorm;
2131 ncoef = INCLENG(ncoef);
2132 }
2133 else if ( *t == NOINDEX ) t++;
2134 else pind[nind++] = *t++;
2135 } while ( t < r );
2136 break;
2137 case SUBEXPRESSION :
2138 if ( t[3] == 0 ) break;
2139 case EXPRESSION :
2140 goto RegEnd;
2141 case ROOTFUNCTION :
2142/*
2143 Tries to take the n-th root inside the rationals
2144 If this is not possible, it clears all flags and
2145 hence tries no more.
2146 Notation:
2147 root_(power(=integer),(rational)number)
2148*/
2149 { WORD nc;
2150 if ( t[2] == 0 ) goto defaultcase;
2151 if ( t[FUNHEAD] != -SNUMBER || t[FUNHEAD+1] < 0 ) goto defaultcase;
2152 if ( t[FUNHEAD+2] == -SNUMBER ) {
2153 if ( t[FUNHEAD+1] == 0 && t[FUNHEAD+3] == 0 ) goto NormZZ;
2154 if ( t[FUNHEAD+1] == 0 ) break;
2155 if ( t[FUNHEAD+3] < 0 ) {
2156 AT.WorkPointer[0] = -t[FUNHEAD+3];
2157 nc = -1;
2158 }
2159 else {
2160 AT.WorkPointer[0] = t[FUNHEAD+3];
2161 nc = 1;
2162 }
2163 AT.WorkPointer[1] = 1;
2164 }
2165 else if ( t[FUNHEAD+2] == t[1]-FUNHEAD-2
2166 && t[FUNHEAD+2] == t[FUNHEAD+2+ARGHEAD]+ARGHEAD
2167 && ABS(t[t[1]-1]) == t[FUNHEAD+2+ARGHEAD] - 1 ) {
2168 WORD *r1, *r2;
2169 if ( t[FUNHEAD+1] == 0 ) break;
2170 i = t[t[1]-1]; r1 = t + FUNHEAD+ARGHEAD+3;
2171 nc = REDLENG(i);
2172 i = ABS(i) - 1;
2173 r2 = AT.WorkPointer;
2174 while ( --i >= 0 ) *r2++ = *r1++;
2175 }
2176 else goto defaultcase;
2177 if ( TakeRatRoot((UWORD *)AT.WorkPointer,&nc,t[FUNHEAD+1]) ) {
2178 t[2] = 0;
2179 goto defaultcase;
2180 }
2181 ncoef = REDLENG(ncoef);
2182 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,nc) )
2183 goto FromNorm;
2184 if ( nc < 0 ) nc = -nc;
2185 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(AT.WorkPointer+nc),nc) )
2186 goto FromNorm;
2187 ncoef = INCLENG(ncoef);
2188 }
2189 break;
2190 case RANDOMFUNCTION :
2191 {
2192 WORD nnc, nc, nca, nr;
2193 UWORD xx;
2194/*
2195 Needs one positive integer argument.
2196 returns (wranf()%argument)+1.
2197 We may call wranf several times to paste UWORDS together
2198 when we need long numbers.
2199 We make little errors when taking the % operator
2200 (not 100% uniform). We correct for that by redoing the
2201 calculation in the (unlikely) case that we are in leftover area
2202*/
2203 if ( t[1] == FUNHEAD ) goto defaultcase;
2204 if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER &&
2205 t[FUNHEAD+1] > 0 ) {
2206 if ( t[FUNHEAD+1] == 1 ) break;
2207redoshort:
2208 ((UWORD *)AT.WorkPointer)[0] = wranf(BHEAD0);
2209 ((UWORD *)AT.WorkPointer)[1] = wranf(BHEAD0);
2210 nr = 2;
2211 if ( ((UWORD *)AT.WorkPointer)[1] == 0 ) {
2212 nr = 1;
2213 if ( ((UWORD *)AT.WorkPointer)[0] == 0 ) {
2214 nr = 0;
2215 }
2216 }
2217 xx = (UWORD)(t[FUNHEAD+1]);
2218 if ( nr ) {
2219 DivLong((UWORD *)AT.WorkPointer,nr
2220 ,&xx,1
2221 ,((UWORD *)AT.WorkPointer)+4,&nnc
2222 ,((UWORD *)AT.WorkPointer)+2,&nc);
2223 ((UWORD *)AT.WorkPointer)[4] = 0;
2224 ((UWORD *)AT.WorkPointer)[5] = 0;
2225 ((UWORD *)AT.WorkPointer)[6] = 1;
2226 DivLong((UWORD *)AT.WorkPointer+4,3
2227 ,&xx,1
2228 ,((UWORD *)AT.WorkPointer)+9,&nnc
2229 ,((UWORD *)AT.WorkPointer)+7,&nca);
2230 AddLong((UWORD *)AT.WorkPointer+4,3
2231 ,((UWORD *)AT.WorkPointer)+7,-nca
2232 ,((UWORD *)AT.WorkPointer)+9,&nnc);
2233 if ( BigLong((UWORD *)AT.WorkPointer,nr
2234 ,((UWORD *)AT.WorkPointer)+9,nnc) >= 0 ) goto redoshort;
2235 }
2236 else nc = 0;
2237 if ( nc == 0 ) {
2238 AT.WorkPointer[2] = (WORD)xx;
2239 nc = 1;
2240 }
2241 ncoef = REDLENG(ncoef);
2242 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+2,nc) )
2243 goto FromNorm;
2244 ncoef = INCLENG(ncoef);
2245 }
2246 else if ( t[FUNHEAD] > 0 && t[1] == t[FUNHEAD]+FUNHEAD
2247 && ABS(t[t[1]-1]) == t[FUNHEAD]-1-ARGHEAD && t[t[1]-1] > 0 ) {
2248 WORD nna, nnb, nni, nnb2, nnb2a;
2249 UWORD *nnt;
2250 nna = t[t[1]-1];
2251 nnb2 = nna-1;
2252 nnb = nnb2/2;
2253 nnt = (UWORD *)(t+t[1]-1-nnb); /* start of denominator */
2254 if ( *nnt != 1 ) goto defaultcase;
2255 for ( nni = 1; nni < nnb; nni++ ) {
2256 if ( nnt[nni] != 0 ) goto defaultcase;
2257 }
2258 nnt = (UWORD *)(t + FUNHEAD + ARGHEAD + 1);
2259
2260 for ( nni = 0; nni < nnb2; nni++ ) {
2261 ((UWORD *)AT.WorkPointer)[nni] = wranf(BHEAD0);
2262 }
2263 nnb2a = nnb2;
2264 while ( nnb2a > 0 && ((UWORD *)AT.WorkPointer)[nnb2a-1] == 0 ) nnb2a--;
2265 if ( nnb2a > 0 ) {
2266 DivLong((UWORD *)AT.WorkPointer,nnb2a
2267 ,nnt,nnb
2268 ,((UWORD *)AT.WorkPointer)+2*nnb2,&nnc
2269 ,((UWORD *)AT.WorkPointer)+nnb2,&nc);
2270 for ( nni = 0; nni < nnb2; nni++ ) {
2271 ((UWORD *)AT.WorkPointer)[nni+2*nnb2] = 0;
2272 }
2273 ((UWORD *)AT.WorkPointer)[3*nnb2] = 1;
2274 DivLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2275 ,nnt,nnb
2276 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc
2277 ,((UWORD *)AT.WorkPointer)+3*nnb2+1,&nca);
2278 AddLong((UWORD *)AT.WorkPointer+2*nnb2,nnb2+1
2279 ,((UWORD *)AT.WorkPointer)+3*nnb2+1,-nca
2280 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,&nnc);
2281 if ( BigLong((UWORD *)AT.WorkPointer,nnb2a
2282 ,((UWORD *)AT.WorkPointer)+4*nnb2+1,nnc) >= 0 ) goto redoshort;
2283 }
2284 else nc = 0;
2285 if ( nc == 0 ) {
2286 for ( nni = 0; nni < nnb; nni++ ) {
2287 ((UWORD *)AT.WorkPointer)[nnb2+nni] = nnt[nni];
2288 }
2289 nc = nnb;
2290 }
2291 ncoef = REDLENG(ncoef);
2292 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,((UWORD *)(AT.WorkPointer))+nnb2,nc) )
2293 goto FromNorm;
2294 ncoef = INCLENG(ncoef);
2295 }
2296 else goto defaultcase;
2297 }
2298 break;
2299 case RANPERM :
2300 if ( *t == RANPERM && t[1] > FUNHEAD && t[FUNHEAD] <= -FUNCTION ) {
2301 WORD **pwork;
2302 WORD *mm, *ww, *ow = AT.WorkPointer;
2303 WORD *Array, *targ, *argstop, narg = 0, itot;
2304 int ie;
2305 argstop = t+t[1];
2306 targ = t+FUNHEAD+1;
2307 while ( targ < argstop ) {
2308 narg++; NEXTARG(targ);
2309 }
2310 WantAddPointers(narg);
2311 pwork = AT.pWorkSpace + AT.pWorkPointer;
2312 targ = t+FUNHEAD+1; narg = 0;
2313 while ( targ < argstop ) {
2314 pwork[narg++] = targ;
2315 NEXTARG(targ);
2316 }
2317/*
2318 Make a random permutation of the numbers 0,...,narg-1
2319 The following code works also for narg == 0 and narg == 1
2320*/
2321 ow = AT.WorkPointer;
2322 Array = AT.WorkPointer;
2323 AT.WorkPointer += narg;
2324 for ( i = 0; i < narg; i++ ) Array[i] = i;
2325 for ( i = 2; i <= narg; i++ ) {
2326 itot = (WORD)(iranf(BHEAD i));
2327 for ( j = 0; j < itot; j++ ) CYCLE1(WORD,Array,i)
2328 }
2329 mm = AT.WorkPointer;
2330 *mm++ = -t[FUNHEAD];
2331 *mm++ = t[1] - 1;
2332 for ( ie = 2; ie < FUNHEAD; ie++ ) *mm++ = t[ie];
2333 for ( i = 0; i < narg; i++ ) {
2334 ww = pwork[Array[i]];
2335 CopyArg(mm,ww);
2336 }
2337 mm = AT.WorkPointer; t++; ww = t;
2338 i = mm[1]; NCOPY(ww,mm,i)
2339 AT.WorkPointer = ow;
2340 goto TryAgain;
2341 }
2342 pnco[nnco++] = t;
2343 break;
2344 case PUTFIRST : /* First argument should be a function, second a number */
2345 if ( ( t[2] & DIRTYFLAG ) != 0 && t[FUNHEAD] <= -FUNCTION
2346 && t[FUNHEAD+1] == -SNUMBER && t[FUNHEAD+2] > 0 ) {
2347 WORD *rr = t+t[1], *mm = t+FUNHEAD+3, *tt, *tt1, *tt2, num = 0;
2348/*
2349 now count the arguments. If not enough: no action.
2350*/
2351 while ( mm < rr ) { num++; NEXTARG(mm); }
2352 if ( num < t[FUNHEAD+2] ) { pnco[nnco++] = t; break; }
2353 /* Replace "putfirst_" with the resulting function code. mm goes to arg start. */
2354 *t = -t[FUNHEAD]; mm = t+FUNHEAD+3;
2355 /* Set i to the arg number we are putting first, then move mm to its start. */
2356 i = t[FUNHEAD+2];
2357 while ( --i > 0 ) { NEXTARG(mm); }
2358 /* Keep a pointer to the arguments trailing the selected: */
2359 WORD *argTail = mm;
2360 NEXTARG(argTail);
2361 tt = TermMalloc("Select_"); /* Move selected out of the way into tmp space */
2362 tt1 = tt;
2363 /* Normal argument: */
2364 if ( *mm > 0 ) {
2365 for ( i = 0; i < *mm; i++ ) *tt1++ = mm[i];
2366 }
2367 /* Fast-notation function: single word */
2368 else if ( *mm <= -FUNCTION ) { *tt1++ = *mm; }
2369 /* Fast-notation symbol, number etc: two words */
2370 else { *tt1++ = mm[0]; *tt1++ = mm[1]; }
2371 /* Put tt2 at the start of the original arguments */
2372 tt2 = t+FUNHEAD+3;
2373 /* Copy leading original arguments after the new first, in the tmp space. */
2374 while ( tt2 < mm ) *tt1++ = *tt2++;
2375 /* i contains the size so far. tt1 goes to the start of the tmp space.
2376 tt2 to the final argument location (we overwrite "putfirst_") */
2377 i = tt1-tt; tt1 = tt; tt2 = t+FUNHEAD;
2378 /* Copy everything so far to its final place */
2379 NCOPY(tt2,tt1,i);
2380 /* We are finished with the tmp space */
2381 TermFree(tt,"Select_");
2382 /* Now copy the trailing args. Use the stored pointer, *mm has been edited
2383 during the copy to tt2 above! NEXTARG(mm) would produce nonsense. */
2384 while ( argTail < rr ) *tt2++ = *argTail++;
2385 /* Set the size of all function args */
2386 t[1] = tt2 - t;
2387 /* Now copy the rest of the term, and set the final term size */
2388 rr = term + *term;
2389 while ( argTail < rr ) *tt2++ = *argTail++;
2390 *term = tt2-term;
2391 if ( functions[*t-FUNCTION].spec == TENSORFUNCTION ) {
2392 // If the output function is a tensor, we have one more job to do.
2393 // The args are formatted as single words, representing indices or vectors.
2394 // There is no ARGHEAD.
2395 WORD *dst = t + FUNHEAD;
2396 WORD *src = dst + 1;
2397 // Strip the type information from the args:
2398 while ( src < t + t[1] ) {
2399 *dst = *src;
2400 dst++; src++; src++;
2401 }
2402 t[1] = dst - t;
2403 // Now copy the rest of the term. We've advanced src one too many above.
2404 src--;
2405 while ( src < term + *term ) {
2406 *dst++ = *src++;
2407 }
2408 *term = dst - term;
2409 }
2410 goto Restart;
2411 }
2412 else pnco[nnco++] = t;
2413 break;
2414#ifdef WITHFLOAT
2415 case FLOATFUN :
2416/*
2417 If it is a proper float_ we give it special treatment.
2418 If it is not proper, we treat it as a regular commuting function.
2419*/
2420 if ( withfloat == 0 ) {
2421 if ( TestFloat(t) == 0 ) goto defaultcase;
2422 firstfloat = t;
2423 withfloat = 1;
2424 }
2425 else { /* There are now at least two float_'s. Time for action. */
2426 if ( TestFloat(t) == 0 ) goto defaultcase;
2427 if ( withfloat == 1 ) UnpackFloat(aux4,firstfloat);
2428 withfloat++;
2429 UnpackFloat(aux5,t);
2430 mpf_mul(aux4,aux4,aux5);
2431 }
2432 break;
2433#endif
2434 case INTFUNCTION :
2435/*
2436 Can be resolved if the first argument is a number
2437 and the second argument either doesn't exist or has
2438 the value +1, 0, -1
2439 +1 : rounding up
2440 0 : rounding towards zero
2441 -1 : rounding down (same as no argument)
2442*/
2443 if ( t[1] <= FUNHEAD ) break;
2444 {
2445 WORD *rr, den, num;
2446 to = t + FUNHEAD;
2447 if ( *to > 0 ) {
2448 if ( *to == ARGHEAD ) goto NormZero;
2449 rr = to + *to;
2450 i = rr[-1];
2451 j = ABS(i);
2452 if ( to[ARGHEAD] != j+1 ) goto NoInteg;
2453 if ( rr >= r ) k = -1;
2454 else if ( *rr == ARGHEAD ) { k = 0; rr += ARGHEAD; }
2455 else if ( *rr == -SNUMBER ) { k = rr[1]; rr += 2; }
2456 else goto NoInteg;
2457 if ( rr != r ) goto NoInteg;
2458 if ( k > 1 || k < -1 ) goto NoInteg;
2459 to += ARGHEAD+1;
2460 j = (j-1) >> 1;
2461 i = ( i < 0 ) ? -j: j;
2462 UnPack((UWORD *)to,i,&den,&num);
2463/*
2464 Potentially the use of NoScrat2 is unsafe.
2465 It makes the routine not reentrant, but because it is
2466 used only locally and because we only call the
2467 low level routines DivLong and AddLong which never
2468 make calls involving Normalize, things are OK after all
2469*/
2470 if ( AN.NoScrat2 == 0 ) {
2471 AN.NoScrat2 = (UWORD *)Malloc1((AM.MaxTal+2)*sizeof(UWORD),"Normalize");
2472 }
2473 if ( DivLong((UWORD *)to,num,(UWORD *)(to+j),den
2474 ,(UWORD *)AT.WorkPointer,&num,AN.NoScrat2,&den) ) goto FromNorm;
2475 if ( k < 0 && den < 0 ) {
2476 *AN.NoScrat2 = 1;
2477 den = -1;
2478 if ( AddLong((UWORD *)AT.WorkPointer,num
2479 ,AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2480 goto FromNorm;
2481 }
2482 else if ( k > 0 && den > 0 ) {
2483 *AN.NoScrat2 = 1;
2484 den = 1;
2485 if ( AddLong((UWORD *)AT.WorkPointer,num,
2486 AN.NoScrat2,den,(UWORD *)AT.WorkPointer,&num) )
2487 goto FromNorm;
2488 }
2489
2490 }
2491 else if ( *to == -SNUMBER ) { /* No rounding needed */
2492 if ( to[1] < 0 ) { *AT.WorkPointer = -to[1]; num = -1; }
2493 else if ( to[1] == 0 ) goto NormZero;
2494 else { *AT.WorkPointer = to[1]; num = 1; }
2495 }
2496 else goto NoInteg;
2497 if ( num == 0 ) goto NormZero;
2498 ncoef = REDLENG(ncoef);
2499 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,num) )
2500 goto FromNorm;
2501 ncoef = INCLENG(ncoef);
2502 break;
2503 }
2504NoInteg:;
2505/*
2506 Fall through if it cannot be resolved
2507*/
2508 default :
2509defaultcase:;
2510 if ( *t < FUNCTION ) {
2511/* INTERNAL_ERROR_EXCL_START */
2512 MLOCK(ErrorMessageLock);
2513 MesPrint("!>Illegal code in Norm");
2514#ifdef DEBUGON
2515 {
2516 UBYTE OutBuf[140];
2517 AO.OutFill = AO.OutputLine = OutBuf;
2518 t = term;
2519 AO.OutSkip = 3;
2520 FiniLine();
2521 i = *t;
2522 while ( --i >= 0 ) {
2523 TalToLine((UWORD)(*t++));
2524 TokenToLine((UBYTE *)" ");
2525 }
2526 AO.OutSkip = 0;
2527 FiniLine();
2528 }
2529#endif
2530 MUNLOCK(ErrorMessageLock);
2531 goto NormMin;
2532/* INTERNAL_ERROR_EXCL_STOP */
2533 }
2534 if ( *t == REPLACEMENT ) {
2535 if ( AR.Eside != LHSIDE ) ReplaceVeto--;
2536 pcom[ncom++] = t;
2537 break;
2538 }
2539/*
2540 if ( *t == AM.termfunnum && t[1] == FUNHEAD+2
2541 && t[FUNHEAD] == -DOLLAREXPRESSION ) termflag++;
2542*/
2543 if ( *t == DUMMYFUN || *t == DUMMYTEN ) {}
2544 else {
2545 if ( *t < (FUNCTION + WILDOFFSET) ) {
2546 if ( ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2547 || ( functions[*t-FUNCTION].minnumargs > 0 ) )
2548 && ( ( t[2] & DIRTYFLAG ) != 0 ) ) {
2549/*
2550 Number of arguments is bounded. And we have not checked.
2551*/
2552 WORD *ta = t + FUNHEAD, *tb = t + t[1];
2553 int numarg = 0;
2554 while ( ta < tb ) { numarg++; NEXTARG(ta) }
2555 if ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2556 && ( numarg >= functions[*t-FUNCTION].maxnumargs ) )
2557 goto NormZero;
2558 if ( ( functions[*t-FUNCTION].minnumargs > 0 )
2559 && ( numarg < functions[*t-FUNCTION].minnumargs ) )
2560 goto NormZero;
2561 }
2562doflags:
2563 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION].tabl == 0 ) ) {
2564 t[2] &= ~DIRTYFLAG;
2565 t[2] |= DIRTYSYMFLAG;
2566 }
2567 if ( functions[*t-FUNCTION].commute ) { pnco[nnco++] = t; }
2568 else { pcom[ncom++] = t; }
2569 }
2570 else {
2571 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION-WILDOFFSET].tabl == 0 ) ) {
2572 t[2] &= ~DIRTYFLAG;
2573 t[2] |= DIRTYSYMFLAG;
2574 }
2575 if ( functions[*t-FUNCTION-WILDOFFSET].commute ) {
2576 pnco[nnco++] = t;
2577 }
2578 else { pcom[ncom++] = t; }
2579 }
2580 }
2581
2582 /* Now hunt for contractible indices */
2583
2584 if ( ( *t < (FUNCTION + WILDOFFSET)
2585 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) || (
2586 *t >= (FUNCTION + WILDOFFSET)
2587 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION ) ) {
2588 if ( *t >= GAMMA && *t <= GAMMASEVEN ) t++;
2589 t += FUNHEAD;
2590 while ( t < r ) {
2591 if ( *t == AM.vectorzero ) goto NormZero;
2592 if ( *t >= AM.OffsetIndex && ( *t >= AM.DumInd
2593 || ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2594 pcon[ncon++] = t;
2595 }
2596 else if ( *t == FUNNYWILD ) { t++; }
2597 t++;
2598 }
2599 }
2600 else {
2601 t += FUNHEAD;
2602 while ( t < r ) {
2603 if ( *t > 0 ) {
2604/*
2605 Here we should worry about a recursion
2606 A problem is the possibility of a construct
2607 like f(mu+nu)
2608*/
2609 t += *t;
2610 }
2611 else if ( *t <= -FUNCTION ) t++;
2612 else if ( *t == -INDEX ) {
2613 if ( t[1] >= AM.OffsetIndex &&
2614 ( t[1] >= AM.DumInd || ( t[1] < AM.WilInd
2615 && indices[t[1]-AM.OffsetIndex].dimension ) ) )
2616 pcon[ncon++] = t+1;
2617 t += 2;
2618 }
2619 else if ( *t == -SYMBOL ) {
2620 if ( t[1] >= MAXPOWER && t[1] < 2*MAXPOWER ) {
2621 *t = -SNUMBER;
2622 t[1] -= MAXPOWER;
2623 }
2624 else if ( t[1] < -MAXPOWER && t[1] > -2*MAXPOWER ) {
2625 *t = -SNUMBER;
2626 t[1] += MAXPOWER;
2627 }
2628 else t += 2;
2629 }
2630 else t += 2;
2631 }
2632 }
2633 break;
2634 }
2635 t = r;
2636TryAgain:;
2637 } while ( t < m );
2638 if ( ANsc ) {
2639 AN.cTerm = ANsc;
2640 r = t = ANsr; m = ANsm;
2641 ANsc = ANsm = ANsr = 0;
2642 goto conscan;
2643 }
2644/*
2645 #] First scan :
2646 #[ Easy denominators :
2647
2648 Easy denominators are denominators that can be replaced by
2649 negative powers of individual subterms. This may add to all
2650 our sublists.
2651
2652*/
2653 if ( nden ) {
2654 for ( k = 0, i = 0; i < nden; i++ ) {
2655 t = pden[i];
2656 if ( ( t[2] & DIRTYFLAG ) == 0 ) continue;
2657 r = t + t[1]; m = t + FUNHEAD;
2658 if ( m >= r ) {
2659 for ( j = i+1; j < nden; j++ ) pden[j-1] = pden[j];
2660 nden--;
2661 for ( j = 0; j < nnco; j++ ) if ( pnco[j] == t ) break;
2662 for ( j++; j < nnco; j++ ) pnco[j-1] = pnco[j];
2663 nnco--;
2664 i--;
2665 }
2666 else {
2667 NEXTARG(m);
2668 if ( m >= r ) continue;
2669/*
2670 We have more than one argument. Split the function.
2671*/
2672 if ( k == 0 ) {
2673 k = 1; to = termout; from = term;
2674 }
2675 while ( from < t ) *to++ = *from++;
2676 m = t + FUNHEAD;
2677 while ( m < r ) {
2678 stop = to;
2679 *to++ = DENOMINATOR;
2680 for ( j = 1; j < FUNHEAD; j++ ) *to++ = 0;
2681 if ( *m < -FUNCTION ) *to++ = *m++;
2682 else if ( *m < 0 ) { *to++ = *m++; *to++ = *m++; }
2683 else {
2684 j = *m; while ( --j >= 0 ) *to++ = *m++;
2685 }
2686 stop[1] = WORDDIF(to,stop);
2687 }
2688 from = r;
2689 if ( i == nden - 1 ) {
2690 stop = term + *term;
2691 while ( from < stop ) *to++ = *from++;
2692 i = *termout = WORDDIF(to,termout);
2693 to = term; from = termout;
2694 while ( --i >= 0 ) *to++ = *from++;
2695 goto Restart;
2696 }
2697 }
2698 }
2699 for ( i = 0; i < nden; i++ ) {
2700 t = pden[i];
2701 if ( ( t[2] & DIRTYFLAG ) == 0 ) continue;
2702 /* Before marking clean, we have to check the denominator argument is
2703 * itself clean and does not contain subexpressions. */
2704 if ( ( t[FUNHEAD] < 0 || t[FUNHEAD+1] == 0 ) && ! TestFunFlag(BHEAD t) ) {
2705 t[2] = 0;
2706 }
2707 if ( t[FUNHEAD] == -SYMBOL ) {
2708 WORD change;
2709 t += FUNHEAD+1;
2710 change = ExtraSymbol(*t,-1,nsym,ppsym,&ncoef);
2711 nsym += change;
2712 ppsym += change * 2;
2713 goto DropDen;
2714 }
2715 else if ( t[FUNHEAD] == -SNUMBER ) {
2716 t += FUNHEAD+1;
2717 if ( *t == 0 ) goto NormInf;
2718 if ( *t < 0 ) { *AT.WorkPointer = -*t; j = -1; }
2719 else { *AT.WorkPointer = *t; j = 1; }
2720 ncoef = REDLENG(ncoef);
2721 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,j) )
2722 goto FromNorm;
2723 ncoef = INCLENG(ncoef);
2724 goto DropDen;
2725 }
2726 else if ( t[FUNHEAD] == ARGHEAD ) goto NormInf;
2727 else if ( t[FUNHEAD] > 0 && t[FUNHEAD+ARGHEAD] ==
2728 t[FUNHEAD]-ARGHEAD ) {
2729 /* Only one term */
2730 r = t + t[1] - 1;
2731 t += FUNHEAD + ARGHEAD + 1;
2732 j = *r;
2733 m = r - ABS(*r) + 1;
2734 if ( j != 3 || ( ( *m != 1 ) || ( m[1] != 1 ) ) ) {
2735 ncoef = REDLENG(ncoef);
2736 if ( DivRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,REDLENG(j),(UWORD *)n_coef,&ncoef) ) goto FromNorm;
2737 ncoef = INCLENG(ncoef);
2738 j = ABS(j) - 3;
2739 t[-FUNHEAD-ARGHEAD] -= j;
2740 t[-ARGHEAD-1] -= j;
2741 t[-1] -= j;
2742 m[0] = m[1] = 1;
2743 m[2] = 3;
2744 }
2745 while ( t < m ) {
2746 r = t + t[1];
2747 if ( *t == SYMBOL || *t == DOTPRODUCT ) {
2748 k = t[1];
2749 pden[i][1] -= k;
2750 pden[i][FUNHEAD] -= k;
2751 pden[i][FUNHEAD+ARGHEAD] -= k;
2752 m -= k;
2753 stop = m + 3;
2754 tt = to = t;
2755 from = r;
2756 if ( *t == SYMBOL ) {
2757 t += 2;
2758 while ( t < r ) {
2759 WORD change;
2760 change = ExtraSymbol(*t,-t[1],nsym,ppsym,&ncoef);
2761 nsym += change;
2762 ppsym += change * 2;
2763 t += 2;
2764 }
2765 }
2766 else {
2767 t += 2;
2768 while ( t < r ) {
2769 *ppdot++ = *t++;
2770 *ppdot++ = *t++;
2771 *ppdot++ = -*t++;
2772 ndot++;
2773 }
2774 }
2775 while ( to < stop ) *to++ = *from++;
2776 r = tt;
2777 }
2778#ifdef WITHFLOAT
2779 else if ( *t == FLOATFUN && TestFloat(t) ) {
2780 k = t[1];
2781 pden[i][1] -= k;
2782 pden[i][FUNHEAD] -= k;
2783 pden[i][FUNHEAD+ARGHEAD] -= k;
2784 UnpackFloat(aux5,t);
2785 if ( withfloat == 0 ) {
2786 mpf_ui_div(aux4,1,aux5);
2787 withfloat = 2;
2788 }
2789 else {
2790 if ( withfloat == 1 ) UnpackFloat(aux4,firstfloat);
2791 mpf_div(aux4,aux4,aux5);
2792 withfloat++;
2793 }
2794 tt = to = t;
2795 while ( r < m ) *to++ = *r++;
2796 *to++ = 1; *to++ = 1; *to++ = 3;
2797 if ( ncoef < 0 ) t[-1] = -t[-1];
2798 m -= k;
2799 stop = m + 3;
2800 r = tt;
2801 }
2802#endif
2803 t = r;
2804 }
2805 if ( pden[i][1] == 4+FUNHEAD+ARGHEAD ) {
2806DropDen:
2807 for ( j = 0; j < nnco; j++ ) {
2808 if ( pden[i] == pnco[j] ) {
2809 --nnco;
2810 while ( j < nnco ) {
2811 pnco[j] = pnco[j+1];
2812 j++;
2813 }
2814 break;
2815 }
2816 }
2817 pden[i--] = pden[--nden];
2818 }
2819 }
2820 }
2821 }
2822/*
2823 #] Easy denominators :
2824 #[ Index Contractions :
2825*/
2826 if ( ndel ) {
2827 t = pdel;
2828 for ( i = 0; i < ndel; i += 2 ) {
2829 if ( t[0] == t[1] ) {
2830 if ( t[0] == EMPTYINDEX ) {}
2831 else if ( *t < AM.OffsetIndex ) {
2832 k = AC.FixIndices[*t];
2833 if ( k < 0 ) { j = -1; k = -k; }
2834 else if ( k > 0 ) j = 1;
2835 else goto NormZero;
2836 goto WithFix;
2837 }
2838 else if ( *t >= AM.DumInd ) {
2839 k = AC.lDefDim;
2840 if ( k ) goto docontract;
2841 }
2842 else if ( *t >= AM.WilInd ) {
2843 k = indices[*t-AM.OffsetIndex-WILDOFFSET].dimension;
2844 if ( k ) goto docontract;
2845 }
2846 else if ( ( k = indices[*t-AM.OffsetIndex].dimension ) != 0 ) {
2847docontract:
2848 if ( k > 0 ) {
2849 j = 1;
2850WithFix: shortnum = k;
2851 ncoef = REDLENG(ncoef);
2852 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),j) )
2853 goto FromNorm;
2854 ncoef = INCLENG(ncoef);
2855 }
2856 else {
2857 WORD change;
2858 change = ExtraSymbol((WORD)(-k),(WORD)1,nsym,ppsym,&ncoef);
2859 nsym += change;
2860 ppsym += change * 2;
2861 }
2862 t[1] = pdel[ndel-1];
2863 t[0] = pdel[ndel-2];
2864HaveCon:
2865 ndel -= 2;
2866 i -= 2;
2867 }
2868 }
2869 else {
2870 if ( *t < AM.OffsetIndex && t[1] < AM.OffsetIndex ) goto NormZero;
2871 j = *t - AM.OffsetIndex;
2872 if ( j >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2873 || ( *t < AM.WilInd && indices[j].dimension ) ) ) {
2874 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2875 if ( *t == *m ) {
2876 *t = m[1];
2877 *m++ = pdel[ndel-2];
2878 *m = pdel[ndel-1];
2879 goto HaveCon;
2880 }
2881 else if ( *t == m[1] ) {
2882 *t = *m;
2883 *m++ = pdel[ndel-2];
2884 *m = pdel[ndel-1];
2885 goto HaveCon;
2886 }
2887 }
2888 }
2889 j = t[1]-AM.OffsetIndex;
2890 if ( j >= 0 && ( ( t[1] >= AM.DumInd && AC.lDefDim )
2891 || ( t[1] < AM.WilInd && indices[j].dimension ) ) ) {
2892 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2893 if ( t[1] == *m ) {
2894 t[1] = m[1];
2895 *m++ = pdel[ndel-2];
2896 *m = pdel[ndel-1];
2897 goto HaveCon;
2898 }
2899 else if ( t[1] == m[1] ) {
2900 t[1] = *m;
2901 *m++ = pdel[ndel-2];
2902 *m = pdel[ndel-1];
2903 goto HaveCon;
2904 }
2905 }
2906 }
2907 t += 2;
2908 }
2909 }
2910 if ( ndel > 0 ) {
2911 if ( nvec ) {
2912 t = pdel;
2913 for ( i = 0; i < ndel; i++ ) {
2914 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2915 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2916 r = pvec + 1;
2917 for ( j = 1; j < nvec; j += 2 ) {
2918 if ( *r == *t ) {
2919 if ( i & 1 ) {
2920 *r = t[-1];
2921 *t-- = pdel[--ndel];
2922 i -= 2;
2923 }
2924 else {
2925 *r = t[1];
2926 t[1] = pdel[--ndel];
2927 i--;
2928 }
2929 *t-- = pdel[--ndel];
2930 break;
2931 }
2932 r += 2;
2933 }
2934 }
2935 t++;
2936 }
2937 }
2938 if ( ndel > 0 && ncon ) {
2939 t = pdel;
2940 for ( i = 0; i < ndel; i++ ) {
2941 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2942 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2943 for ( j = 0; j < ncon; j++ ) {
2944 if ( *pcon[j] == *t ) {
2945 if ( i & 1 ) {
2946 *pcon[j] = t[-1];
2947 *t-- = pdel[--ndel];
2948 i -= 2;
2949 }
2950 else {
2951 *pcon[j] = t[1];
2952 t[1] = pdel[--ndel];
2953 i--;
2954 }
2955 *t-- = pdel[--ndel];
2956 didcontr++;
2957 r = pcon[j];
2958 for ( j = 0; j < nnco; j++ ) {
2959 m = pnco[j];
2960 if ( r > m && r < m+m[1] ) {
2961 m[2] |= DIRTYSYMFLAG;
2962 break;
2963 }
2964 }
2965 for ( j = 0; j < ncom; j++ ) {
2966 m = pcom[j];
2967 if ( r > m && r < m+m[1] ) {
2968 m[2] |= DIRTYSYMFLAG;
2969 break;
2970 }
2971 }
2972 for ( j = 0; j < neps; j++ ) {
2973 m = peps[j];
2974 if ( r > m && r < m+m[1] ) {
2975 m[2] |= DIRTYSYMFLAG;
2976 break;
2977 }
2978 }
2979 break;
2980 }
2981 }
2982 }
2983 t++;
2984 }
2985 }
2986 }
2987 }
2988 if ( nvec ) {
2989 t = pvec + 1;
2990 for ( i = 3; i < nvec; i += 2 ) {
2991 k = *t - AM.OffsetIndex;
2992 if ( k >= 0 && ( ( *t > AM.DumInd && AC.lDefDim )
2993 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2994 r = t + 2;
2995 for ( j = i; j < nvec; j += 2 ) {
2996 if ( *r == *t ) { /* Another dotproduct */
2997 *ppdot++ = t[-1];
2998 *ppdot++ = r[-1];
2999 *ppdot++ = 1;
3000 ndot++;
3001 *r-- = pvec[--nvec];
3002 *r = pvec[--nvec];
3003 *t-- = pvec[--nvec];
3004 *t-- = pvec[--nvec];
3005 i -= 2;
3006 break;
3007 }
3008 r += 2;
3009 }
3010 }
3011 t += 2;
3012 }
3013 if ( nvec > 0 && ncon ) {
3014 t = pvec + 1;
3015 for ( i = 1; i < nvec; i += 2 ) {
3016 k = *t - AM.OffsetIndex;
3017 if ( k >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
3018 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
3019 for ( j = 0; j < ncon; j++ ) {
3020 if ( *pcon[j] == *t ) {
3021 *pcon[j] = t[-1];
3022 *t-- = pvec[--nvec];
3023 *t-- = pvec[--nvec];
3024 r = pcon[j];
3025 pcon[j] = pcon[--ncon];
3026 i -= 2;
3027 for ( j = 0; j < nnco; j++ ) {
3028 m = pnco[j];
3029 if ( r > m && r < m+m[1] ) {
3030 m[2] |= DIRTYSYMFLAG;
3031 break;
3032 }
3033 }
3034 for ( j = 0; j < ncom; j++ ) {
3035 m = pcom[j];
3036 if ( r > m && r < m+m[1] ) {
3037 m[2] |= DIRTYSYMFLAG;
3038 break;
3039 }
3040 }
3041 for ( j = 0; j < neps; j++ ) {
3042 m = peps[j];
3043 if ( r > m && r < m+m[1] ) {
3044 m[2] |= DIRTYSYMFLAG;
3045 break;
3046 }
3047 }
3048 break;
3049 }
3050 }
3051 }
3052 t += 2;
3053 }
3054 }
3055 }
3056/*
3057 #] Index Contractions :
3058 #[ NonCommuting Functions :
3059*/
3060 m = fillsetexp;
3061 if ( nnco ) {
3062 for ( i = 0; i < nnco; i++ ) {
3063 t = pnco[i];
3064 if ( ( *t >= (FUNCTION+WILDOFFSET)
3065 && functions[*t-FUNCTION-WILDOFFSET].spec <= 0 )
3066 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3067 && functions[*t-FUNCTION].spec <= 0 ) ) {
3068 DoRevert(t,m);
3069 if ( didcontr ) {
3070 r = t + FUNHEAD;
3071 t += t[1];
3072 while ( r < t ) {
3073 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3074 *r = -SNUMBER;
3075 didcontr--;
3076 pnco[i][2] |= DIRTYSYMFLAG;
3077 }
3078 NEXTARG(r)
3079 }
3080 }
3081 }
3082 }
3083
3084 /* First should come the code for function properties. */
3085
3086 /* First we test for symmetric properties and the DIRTYSYMFLAG */
3087
3088 for ( i = 0; i < nnco; i++ ) {
3089 t = pnco[i];
3090 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) && *t != DOLLAREXPRESSION ) {
3091 l = 0; /* to make the compiler happy */
3092 if ( ( *t >= (FUNCTION+WILDOFFSET)
3093 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3094 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3095 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3096 if ( *t >= (FUNCTION+WILDOFFSET) ) {
3097 *t -= WILDOFFSET;
3098 j = FullSymmetrize(BHEAD t,l);
3099 *t += WILDOFFSET;
3100 }
3101 else j = FullSymmetrize(BHEAD t,l);
3102 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3103 if ( ( j & 2 ) != 0 ) goto NormZero;
3104 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3105 }
3106 }
3107 else t[2] &= ~DIRTYSYMFLAG;
3108 }
3109 }
3110
3111 /* Non commuting functions are then tested for commutation
3112 rules. If needed their order is exchanged. */
3113
3114 k = nnco - 1;
3115 for ( i = 0; i < k; i++ ) {
3116 j = i;
3117 while ( Commute(pnco[j],pnco[j+1]) ) {
3118 t = pnco[j]; pnco[j] = pnco[j+1]; pnco[j+1] = t;
3119 l = j-1;
3120 while ( l >= 0 && Commute(pnco[l],pnco[l+1]) ) {
3121 t = pnco[l]; pnco[l] = pnco[l+1]; pnco[l+1] = t;
3122 l--;
3123 }
3124 if ( ++j >= k ) break;
3125 }
3126 }
3127
3128 /* Finally they are written to output. gamma matrices
3129 are bundled if possible */
3130
3131 for ( i = 0; i < nnco; i++ ) {
3132 t = pnco[i];
3133 if ( *t == IDFUNCTION ) AN.idfunctionflag = 1;
3134 if ( *t >= GAMMA && *t <= GAMMASEVEN ) {
3135 WORD gtype;
3136 to = m;
3137 *m++ = GAMMA;
3138 m++;
3139 FILLFUN(m)
3140 *m++ = stype = t[FUNHEAD]; /* type of string */
3141 j = 0;
3142 nnum = 0;
3143 do {
3144 r = t + t[1];
3145 if ( *t == GAMMAFIVE ) {
3146 gtype = GAMMA5; t += FUNHEAD; goto onegammamatrix; }
3147 else if ( *t == GAMMASIX ) {
3148 gtype = GAMMA6; t += FUNHEAD; goto onegammamatrix; }
3149 else if ( *t == GAMMASEVEN ) {
3150 gtype = GAMMA7; t += FUNHEAD; goto onegammamatrix; }
3151 t += FUNHEAD+1;
3152 while ( t < r ) {
3153 gtype = *t;
3154onegammamatrix:
3155 if ( gtype == GAMMA5 ) {
3156 if ( j == GAMMA1 ) j = GAMMA5;
3157 else if ( j == GAMMA5 ) j = GAMMA1;
3158 else if ( j == GAMMA7 ) ncoef = -ncoef;
3159 if ( nnum & 1 ) ncoef = -ncoef;
3160 }
3161 else if ( gtype == GAMMA6 || gtype == GAMMA7 ) {
3162 if ( nnum & 1 ) {
3163 if ( gtype == GAMMA6 ) gtype = GAMMA7;
3164 else gtype = GAMMA6;
3165 }
3166 if ( j == GAMMA1 ) j = gtype;
3167 else if ( j == GAMMA5 ) {
3168 j = gtype;
3169 if ( j == GAMMA7 ) ncoef = -ncoef;
3170 }
3171 else if ( j != gtype ) goto NormZero;
3172 else {
3173 shortnum = 2;
3174 ncoef = REDLENG(ncoef);
3175 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),1) ) goto FromNorm;
3176 ncoef = INCLENG(ncoef);
3177 }
3178 }
3179 else {
3180 *m++ = gtype; nnum++;
3181 }
3182 t++;
3183 }
3184
3185 } while ( ( ++i < nnco ) && ( *(t = pnco[i]) >= GAMMA
3186 && *t <= GAMMASEVEN ) && ( t[FUNHEAD] == stype ) );
3187 i--;
3188 if ( j ) {
3189 k = WORDDIF(m,to) - FUNHEAD-1;
3190 r = m;
3191 from = m++;
3192 while ( --k >= 0 ) *from-- = *--r;
3193 *from = j;
3194 }
3195 to[1] = WORDDIF(m,to);
3196 }
3197 else if ( *t < 0 ) {
3198 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3199 FILLFUN3(m)
3200 }
3201 else {
3202 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3203 && *t != REPLACEMENT && *t != DOLLAREXPRESSION
3204 && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3205 k = t[1];
3206 NCOPY(m,t,k);
3207 }
3208 }
3209
3210 }
3211/*
3212 #] NonCommuting Functions :
3213 #[ Commuting Functions :
3214*/
3215 if ( ncom ) {
3216 for ( i = 0; i < ncom; i++ ) {
3217 t = pcom[i];
3218 if ( ( *t >= (FUNCTION+WILDOFFSET)
3219 && functions[*t-FUNCTION-WILDOFFSET].spec <= 0 )
3220 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3221 && functions[*t-FUNCTION].spec <= 0 ) ) {
3222 DoRevert(t,m);
3223 if ( didcontr ) {
3224 r = t + FUNHEAD;
3225 t += t[1];
3226 while ( r < t ) {
3227 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3228 *r = -SNUMBER;
3229 didcontr--;
3230 pcom[i][2] |= DIRTYSYMFLAG;
3231 }
3232 NEXTARG(r)
3233 }
3234 }
3235 }
3236 }
3237
3238 /* Now we test for symmetric properties and the DIRTYSYMFLAG */
3239
3240 for ( i = 0; i < ncom; i++ ) {
3241 t = pcom[i];
3242 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) ) {
3243 l = 0; /* to make the compiler happy */
3244 if ( ( *t >= (FUNCTION+WILDOFFSET)
3245 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3246 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3247 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3248 if ( *t >= (FUNCTION+WILDOFFSET) ) {
3249 *t -= WILDOFFSET;
3250 j = FullSymmetrize(BHEAD t,l);
3251 *t += WILDOFFSET;
3252 }
3253 else j = FullSymmetrize(BHEAD t,l);
3254 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3255 if ( ( j & 2 ) != 0 ) goto NormZero;
3256 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3257 }
3258 }
3259 else t[2] &= ~DIRTYSYMFLAG;
3260 }
3261 }
3262/*
3263 Sort the functions
3264 From a purists point of view this can be improved.
3265 There are slow and fast arguments and no conversions are
3266 taken into account here.
3267*/
3268 for ( i = 1; i < ncom; i++ ) {
3269 for ( j = i; j > 0; j-- ) {
3270 WORD jj,kk;
3271 jj = j-1;
3272 t = pcom[jj];
3273 r = pcom[j];
3274 if ( *t < 0 ) {
3275 if ( *r < 0 ) { if ( *t >= *r ) goto NextI; }
3276 else { if ( -*t <= *r ) goto NextI; }
3277 goto jexch;
3278 }
3279 else if ( *r < 0 ) {
3280 if ( *t < -*r ) goto NextI;
3281 goto jexch;
3282 }
3283 else if ( *t != *r ) {
3284 if ( *t < *r ) goto NextI;
3285jexch: t = pcom[j]; pcom[j] = pcom[jj]; pcom[jj] = t;
3286 continue;
3287 }
3288 if ( AC.properorderflag ) {
3289 if ( ( *t >= (FUNCTION+WILDOFFSET)
3290 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
3291 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3292 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) ) {}
3293 else {
3294 WORD *s1, *s2, *ss1, *ss2;
3295 s1 = t+FUNHEAD; s2 = r+FUNHEAD;
3296 ss1 = t + t[1]; ss2 = r + r[1];
3297 while ( s1 < ss1 && s2 < ss2 ) {
3298 k = CompArg(s1,s2);
3299 if ( k > 0 ) goto jexch;
3300 if ( k < 0 ) goto NextI;
3301 NEXTARG(s1)
3302 NEXTARG(s2)
3303 }
3304 if ( s1 < ss1 ) goto jexch;
3305 goto NextI;
3306 }
3307 k = t[1] - FUNHEAD;
3308 kk = r[1] - FUNHEAD;
3309 t += FUNHEAD;
3310 r += FUNHEAD;
3311 while ( k > 0 && kk > 0 ) {
3312 if ( *t < *r ) goto NextI;
3313 else if ( *t++ > *r++ ) goto jexch;
3314 k--; kk--;
3315 }
3316 if ( k > 0 ) goto jexch;
3317 goto NextI;
3318 }
3319 else
3320 {
3321 k = t[1] - FUNHEAD;
3322 kk = r[1] - FUNHEAD;
3323 t += FUNHEAD;
3324 r += FUNHEAD;
3325 while ( k > 0 && kk > 0 ) {
3326 if ( *t < *r ) goto NextI;
3327 else if ( *t++ > *r++ ) goto jexch;
3328 k--; kk--;
3329 }
3330 if ( k > 0 ) goto jexch;
3331 goto NextI;
3332 }
3333 }
3334NextI:;
3335 }
3336 for ( i = 0; i < ncom; i++ ) {
3337 t = pcom[i];
3338 if ( *t == THETA || *t == THETA2 ) {
3339 if ( ( k = DoTheta(BHEAD t) ) == 0 ) goto NormZero;
3340 else if ( k < 0 ) {
3341 k = t[1];
3342 NCOPY(m,t,k);
3343 }
3344 }
3345 else if ( *t == DELTA2 || *t == DELTAP ) {
3346 if ( ( k = DoDelta(t) ) == 0 ) goto NormZero;
3347 else if ( k < 0 ) {
3348 k = t[1];
3349 NCOPY(m,t,k);
3350 }
3351 }
3352 else if ( *t == AR.PolyFunInv && AR.PolyFunType == 2 ) {
3353/*
3354 If there are two arguments, exchange them, change the
3355 name of the function and go to dealing with PolyRatFun.
3356*/
3357 WORD *mm, *tt = t, numt = 0;
3358 tt += FUNHEAD;
3359 while ( tt < t+t[1] ) { numt++; NEXTARG(tt) }
3360 if ( numt == 2 ) {
3361 tt = t; mm = m; k = t[1];
3362 NCOPY(mm,tt,k)
3363 mm = m+FUNHEAD;
3364 NEXTARG(mm);
3365 tt = t+FUNHEAD;
3366 if ( *mm < 0 ) {
3367 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3368 else { *tt++ = *mm++; *tt++ = *mm++; }
3369 }
3370 else {
3371 k = *mm; NCOPY(tt,mm,k)
3372 }
3373 mm = m+FUNHEAD;
3374 if ( *mm < 0 ) {
3375 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3376 else { *tt++ = *mm++; *tt++ = *mm++; }
3377 }
3378 else {
3379 k = *mm; NCOPY(tt,mm,k)
3380 }
3381 *t = AR.PolyFun;
3382 t[2] |= MUSTCLEANPRF;
3383 goto regularratfun;
3384 }
3385 }
3386 else if ( *t == AR.PolyFun ) {
3387 if ( AR.PolyFunType == 1 ) { /* Regular PolyFun with one argument */
3388 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3389 t[1] == FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) goto NormZero;
3390 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) {
3391 if ( AN.PolyNormFlag == 0 ) {
3392 AN.PolyNormFlag = 1;
3393 AN.PolyFunTodo = 0;
3394 }
3395 }
3396 k = t[1];
3397 NCOPY(m,t,k);
3398 }
3399 else if ( AR.PolyFunType == 2 ) {
3400/*
3401 PolyRatFun.
3402 Regular type: Two arguments
3403 Power expanded: One argument. Here to be treated as
3404 AR.PolyFunType == 1, but with power cutoff.
3405*/
3406regularratfun:;
3407/*
3408 First check for zeroes.
3409*/
3410 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3411 t[1] > FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) {
3412 u = t + FUNHEAD + 2;
3413 if ( *u < 0 ) {
3414 if ( *u <= -FUNCTION ) {}
3415 else if ( t[1] == FUNHEAD+4 && t[FUNHEAD+2] == -SNUMBER
3416 && t[FUNHEAD+3] == 0 ) goto NormPRF;
3417 else if ( t[1] == FUNHEAD+4 ) goto NormZero;
3418 }
3419 else if ( t[1] == *u+FUNHEAD+2 ) goto NormZero;
3420 }
3421 else {
3422 u = t+FUNHEAD; NEXTARG(u);
3423 if ( *u == -SNUMBER && u[1] == 0 ) goto NormInf;
3424 }
3425 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3426 else if ( i < ncom-1 && pcom[i+1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3427 k = t[1];
3428 if ( AN.PolyNormFlag ) {
3429 if ( AR.PolyFunExp == 0 ) {
3430 AN.PolyFunTodo = 0;
3431 NCOPY(m,t,k);
3432 }
3433 else if ( AR.PolyFunExp == 1 ) { /* get highest divergence */
3434 if ( PolyFunMode == 0 ) {
3435 NCOPY(m,t,k);
3436 AN.PolyFunTodo = 1;
3437 }
3438 else {
3439 WORD *mmm = m;
3440 NCOPY(m,t,k);
3441 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3442 goto FromNorm;
3443 m = mmm+mmm[1];
3444 }
3445 }
3446 else {
3447 if ( PolyFunMode == 0 ) {
3448 NCOPY(m,t,k);
3449 AN.PolyFunTodo = 1;
3450 }
3451 else {
3452 WORD *mmm = m;
3453 NCOPY(m,t,k);
3454 if ( ExpandRat(BHEAD mmm) != 0 )
3455 goto FromNorm;
3456 m = mmm+mmm[1];
3457 }
3458 }
3459 }
3460 else {
3461 if ( AR.PolyFunExp == 0 ) {
3462 AN.PolyFunTodo = 0;
3463 NCOPY(m,t,k);
3464 }
3465 else if ( AR.PolyFunExp == 1 ) { /* get highest divergence */
3466 WORD *mmm = m;
3467 NCOPY(m,t,k);
3468 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3469 goto FromNorm;
3470 m = mmm+mmm[1];
3471 }
3472 else {
3473 WORD *mmm = m;
3474 NCOPY(m,t,k);
3475 if ( ExpandRat(BHEAD mmm) != 0 )
3476 goto FromNorm;
3477 m = mmm+mmm[1];
3478 }
3479 }
3480 }
3481 }
3482 else if ( *t > 0 ) {
3483 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3484 && *t != REPLACEMENT && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3485 k = t[1];
3486 NCOPY(m,t,k);
3487 }
3488 else {
3489 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3490 FILLFUN3(m)
3491 }
3492 }
3493 }
3494/*
3495 #] Commuting Functions :
3496 #[ Track Replace_ :
3497*/
3498 if ( ReplaceVeto < 0 ) {
3499/*
3500 We found one (or more) replace_ functions and all other
3501 functions are 'clean' (no dirty flag).
3502 Now we check whether one of these functions can be used.
3503 Thus far the functions go from fillsetexp to m.
3504 Somewhere in there there are -ReplaceVeto occurrences of REPLACEMENT.
3505 Hunt for the first one that fits the bill.
3506 Note that replace_ is a commuting function.
3507*/
3508 WORD *ma = fillsetexp, *mb, *mc;
3509 while ( ma < m ) {
3510 mb = ma + ma[1];
3511 if ( *ma != REPLACEMENT ) {
3512 ma = mb;
3513 continue;
3514 }
3515 if ( *ma == REPLACEMENT && ReplaceType == -1 ) {
3516 mc = ma;
3517 ReplaceType = 0;
3518 if ( AN.RSsize < 2*ma[1]+SUBEXPSIZE ) {
3519 if ( AN.ReplaceScrat ) M_free(AN.ReplaceScrat,"AN.ReplaceScrat");
3520 AN.RSsize = 2*ma[1]+SUBEXPSIZE+40;
3521 AN.ReplaceScrat = (WORD *)Malloc1((AN.RSsize+1)*sizeof(WORD),"AN.ReplaceScrat");
3522 }
3523 ma += FUNHEAD;
3524 ReplaceSub = AN.ReplaceScrat;
3525 ReplaceSub += SUBEXPSIZE;
3526 while ( ma < mb ) {
3527 if ( *ma > 0 ) goto NoRep;
3528 if ( *ma <= -FUNCTION ) {
3529 *ReplaceSub++ = FUNTOFUN;
3530 *ReplaceSub++ = 4;
3531 *ReplaceSub++ = -*ma++;
3532 if ( *ma > -FUNCTION ) goto NoRep;
3533 *ReplaceSub++ = -*ma++;
3534 }
3535 else if ( ma+4 > mb ) goto NoRep;
3536 else {
3537 if ( *ma == -SYMBOL ) {
3538 if ( ma[2] == -SYMBOL && ma+4 <= mb )
3539 *ReplaceSub++ = SYMTOSYM;
3540 else if ( ma[2] == -SNUMBER && ma+4 <= mb ) {
3541 *ReplaceSub++ = SYMTONUM;
3542 if ( ReplaceType == 0 ) {
3543 oldtoprhs = C->numrhs;
3544 oldcpointer = C->Pointer - C->Buffer;
3545 }
3546 ReplaceType = 1;
3547 }
3548 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3549 *ReplaceSub++ = SYMTONUM;
3550 *ReplaceSub++ = 4;
3551 *ReplaceSub++ = ma[1];
3552 *ReplaceSub++ = 0;
3553 ma += 2+ARGHEAD;
3554 continue;
3555 }
3556/*
3557 Next is the subexpression. We have to test that
3558 it isn't vector-like or index-like
3559*/
3560 else if ( ma[2] > 0 ) {
3561 WORD *sstop, *ttstop, n;
3562 ss = ma+2;
3563 sstop = ss + *ss;
3564 ss += ARGHEAD;
3565 while ( ss < sstop ) {
3566 tt = ss + *ss;
3567 ttstop = tt - ABS(tt[-1]);
3568 ss++;
3569 while ( ss < ttstop ) {
3570 if ( *ss == INDEX ) goto NoRep;
3571 ss += ss[1];
3572 }
3573 ss = tt;
3574 }
3575 subtype = SYMTOSUB;
3576 if ( ReplaceType == 0 ) {
3577 oldtoprhs = C->numrhs;
3578 oldcpointer = C->Pointer - C->Buffer;
3579 }
3580 ReplaceType = 1;
3581 ss = AddRHS(AT.ebufnum,1);
3582 tt = ma+2;
3583 n = *tt - ARGHEAD;
3584 tt += ARGHEAD;
3585 while ( (ss + n + 10) > C->Top ) ss = DoubleCbuffer(AT.ebufnum,ss,14);
3586 while ( --n >= 0 ) *ss++ = *tt++;
3587 *ss++ = 0;
3588 C->rhs[C->numrhs+1] = ss;
3589 C->Pointer = ss;
3590 *ReplaceSub++ = subtype;
3591 *ReplaceSub++ = 4;
3592 *ReplaceSub++ = ma[1];
3593 *ReplaceSub++ = C->numrhs;
3594 ma += 2 + ma[2];
3595 continue;
3596 }
3597 else goto NoRep;
3598 }
3599 else if ( ( *ma == -VECTOR || *ma == -MINVECTOR ) && ma+4 <= mb ) {
3600 if ( ma[2] == -VECTOR ) {
3601 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOVEC;
3602 else *ReplaceSub++ = VECTOMIN;
3603 }
3604 else if ( ma[2] == -MINVECTOR ) {
3605 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOMIN;
3606 else *ReplaceSub++ = VECTOVEC;
3607 }
3608/*
3609 Next is a vector-like subexpression
3610 Search for vector nature first
3611*/
3612 else if ( ma[2] > 0 ) {
3613 WORD *sstop, *ttstop, *w, *mm, n, count;
3614 WORD *v1, *v2 = 0;
3615 if ( *ma == -MINVECTOR ) {
3616 ss = ma+2;
3617 sstop = ss + *ss;
3618 ss += ARGHEAD;
3619 while ( ss < sstop ) {
3620 ss += *ss;
3621 ss[-1] = -ss[-1];
3622 }
3623 *ma = -VECTOR;
3624 }
3625 ss = ma+2;
3626 sstop = ss + *ss;
3627 ss += ARGHEAD;
3628 while ( ss < sstop ) {
3629 tt = ss + *ss;
3630 ttstop = tt - ABS(tt[-1]);
3631 ss++;
3632 count = 0;
3633 while ( ss < ttstop ) {
3634 if ( *ss == INDEX ) {
3635 n = ss[1] - 2; ss += 2;
3636 while ( --n >= 0 ) {
3637 if ( *ss < MINSPEC ) count++;
3638 ss++;
3639 }
3640 }
3641 else ss += ss[1];
3642 }
3643 if ( count != 1 ) goto NoRep;
3644 ss = tt;
3645 }
3646 subtype = VECTOSUB;
3647 if ( ReplaceType == 0 ) {
3648 oldtoprhs = C->numrhs;
3649 oldcpointer = C->Pointer - C->Buffer;
3650 }
3651 ReplaceType = 1;
3652 mm = AddRHS(AT.ebufnum,1);
3653 *ReplaceSub++ = subtype;
3654 *ReplaceSub++ = 4;
3655 *ReplaceSub++ = ma[1];
3656 *ReplaceSub++ = C->numrhs;
3657 w = ma+2;
3658 n = *w - ARGHEAD;
3659 w += ARGHEAD;
3660 while ( (mm + n + 10) > C->Top )
3661 mm = DoubleCbuffer(AT.ebufnum,mm,15);
3662 while ( --n >= 0 ) *mm++ = *w++;
3663 *mm++ = 0;
3664 C->rhs[C->numrhs+1] = mm;
3665 C->Pointer = mm;
3666 mm = AddRHS(AT.ebufnum,1);
3667 w = ma+2;
3668 n = *w - ARGHEAD;
3669 w += ARGHEAD;
3670 while ( (mm + n + 13) > C->Top )
3671 mm = DoubleCbuffer(AT.ebufnum,mm,16);
3672 sstop = w + n;
3673 while ( w < sstop ) {
3674 tt = w + *w; ttstop = tt - ABS(tt[-1]);
3675 ss = mm; mm++; w++;
3676 while ( w < ttstop ) { /* Subterms */
3677 if ( *w != INDEX ) {
3678 n = w[1];
3679 NCOPY(mm,w,n);
3680 }
3681 else {
3682 v1 = mm;
3683 *mm++ = *w++;
3684 *mm++ = n = *w++;
3685 n -= 2;
3686 while ( --n >= 0 ) {
3687 if ( *w >= MINSPEC ) *mm++ = *w++;
3688 else v2 = w++;
3689 }
3690 n = WORDDIF(mm,v1);
3691 if ( n != v1[1] ) {
3692 if ( n <= 2 ) mm -= 2;
3693 else v1[1] = n;
3694 *mm++ = VECTOR;
3695 *mm++ = 4;
3696 *mm++ = *v2;
3697 *mm++ = FUNNYVEC;
3698 }
3699 }
3700 }
3701 while ( w < tt ) *mm++ = *w++;
3702 *ss = WORDDIF(mm,ss);
3703 }
3704 *mm++ = 0;
3705 C->rhs[C->numrhs+1] = mm;
3706 C->Pointer = mm;
3707 if ( mm > C->Top ) {
3708/* INTERNAL_ERROR_EXCL_START */
3709 MLOCK(ErrorMessageLock);
3710 MesPrint("!>Internal error in Normalize with extra compiler buffer");
3711 MUNLOCK(ErrorMessageLock);
3712 Terminate(-1);
3713/* INTERNAL_ERROR_EXCL_STOP */
3714 }
3715 ma += 2 + ma[2];
3716 continue;
3717 }
3718 else goto NoRep;
3719 }
3720 else if ( *ma == -INDEX ) {
3721 if ( ( ma[2] == -INDEX || ma[2] == -VECTOR )
3722 && ma+4 <= mb )
3723 *ReplaceSub++ = INDTOIND;
3724 else if ( ma[1] >= AM.OffsetIndex ) {
3725 if ( ma[2] == -SNUMBER && ma+4 <= mb
3726 && ma[3] >= 0 && ma[3] < AM.OffsetIndex )
3727 *ReplaceSub++ = INDTOIND;
3728 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3729 *ReplaceSub++ = INDTOIND;
3730 *ReplaceSub++ = 4;
3731 *ReplaceSub++ = ma[1];
3732 *ReplaceSub++ = 0;
3733 ma += 2+ARGHEAD;
3734 continue;
3735 }
3736 else goto NoRep;
3737 }
3738 else goto NoRep;
3739 }
3740 else goto NoRep;
3741 *ReplaceSub++ = 4;
3742 *ReplaceSub++ = ma[1];
3743 *ReplaceSub++ = ma[3];
3744 ma += 4;
3745 }
3746
3747 }
3748 AN.ReplaceScrat[1] = ReplaceSub-AN.ReplaceScrat;
3749/*
3750 Success. This means that we have to remove the replace_
3751 from the functions. It starts at mc and end at mb.
3752*/
3753 while ( mb < m ) *mc++ = *mb++;
3754 m = mc;
3755 break;
3756NoRep:
3757 if ( ReplaceType > 0 ) {
3758 C->numrhs = oldtoprhs;
3759 C->Pointer = C->Buffer + oldcpointer;
3760 }
3761 ReplaceType = -1;
3762 if ( ++ReplaceVeto >= 0 ) break;
3763 }
3764 ma = mb;
3765 }
3766 }
3767/*
3768 #] Track Replace_ :
3769 #[ LeviCivita tensors :
3770*/
3771 if ( neps ) {
3772 to = m;
3773 for ( i = 0; i < neps; i++ ) { /* Put the indices in order */
3774 t = peps[i];
3775 if ( ( t[2] & DIRTYSYMFLAG ) != DIRTYSYMFLAG ) continue;
3776 t[2] &= ~DIRTYSYMFLAG;
3777 if ( AR.Eside == LHSIDE || AR.Eside == LHSIDEX ) {
3778 /* Potential problems with FUNNYWILD */
3779/*
3780 First make sure all FUNNIES are at the end.
3781 Then sort separately
3782*/
3783 r = t + FUNHEAD;
3784 m = tt = t + t[1];
3785 while ( r < m ) {
3786 if ( *r != FUNNYWILD ) { r++; continue; }
3787 k = r[1]; u = r + 2;
3788 while ( u < tt ) {
3789 u[-2] = *u;
3790 if ( *u != FUNNYWILD ) ncoef = -ncoef;
3791 u++;
3792 }
3793 tt[-2] = FUNNYWILD; tt[-1] = k; m -= 2;
3794 }
3795 t += FUNHEAD;
3796 do {
3797 for ( r = t + 1; r < m; r++ ) {
3798 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3799 else if ( *r == *t ) goto NormZero;
3800 }
3801 t++;
3802 } while ( t < m );
3803 do {
3804 for ( r = t + 2; r < tt; r += 2 ) {
3805 if ( r[1] < t[1] ) {
3806 k = r[1]; r[1] = t[1]; t[1] = k; ncoef = -ncoef; }
3807 else if ( r[1] == t[1] ) goto NormZero;
3808 }
3809 t += 2;
3810 } while ( t < tt );
3811 }
3812 else {
3813 m = t + t[1];
3814 t += FUNHEAD;
3815 do {
3816 for ( r = t + 1; r < m; r++ ) {
3817 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3818 else if ( *r == *t ) goto NormZero;
3819 }
3820 t++;
3821 } while ( t < m );
3822 }
3823 }
3824
3825 /* Sort the tensors */
3826
3827 for ( i = 0; i < (neps-1); i++ ) {
3828 t = peps[i];
3829 for ( j = i+1; j < neps; j++ ) {
3830 r = peps[j];
3831 if ( t[1] > r[1] ) {
3832 peps[i] = m = r; peps[j] = r = t; t = m;
3833 }
3834 else if ( t[1] == r[1] ) {
3835 k = t[1] - FUNHEAD;
3836 m = t + FUNHEAD;
3837 r += FUNHEAD;
3838 do {
3839 if ( *r < *m ) {
3840 m = peps[j]; peps[j] = t; peps[i] = t = m;
3841 break;
3842 }
3843 else if ( *r++ > *m++ ) break;
3844 } while ( --k > 0 );
3845 }
3846 }
3847 }
3848 m = to;
3849 for ( i = 0; i < neps; i++ ) {
3850 t = peps[i];
3851 k = t[1];
3852 NCOPY(m,t,k);
3853 }
3854 }
3855/*
3856 #] LeviCivita tensors :
3857 #[ Delta :
3858*/
3859 if ( ndel ) {
3860 r = t = pdel;
3861 for ( i = 0; i < ndel; i += 2, r += 2 ) {
3862 if ( r[1] < r[0] ) { k = *r; *r = r[1]; r[1] = k; }
3863 }
3864 for ( i = 2; i < ndel; i += 2, t += 2 ) {
3865 r = t + 2;
3866 for ( j = i; j < ndel; j += 2 ) {
3867 if ( *r > *t ) { r += 2; }
3868 else if ( *r < *t ) {
3869 k = *r; *r++ = *t; *t++ = k;
3870 k = *r; *r++ = *t; *t-- = k;
3871 }
3872 else {
3873 if ( *++r < t[1] ) {
3874 k = *r; *r = t[1]; t[1] = k;
3875 }
3876 r++;
3877 }
3878 }
3879 }
3880 t = pdel;
3881 *m++ = DELTA;
3882 *m++ = ndel + 2;
3883 i = ndel;
3884 NCOPY(m,t,i);
3885 }
3886/*
3887 #] Delta :
3888 #[ Loose Vectors/Indices :
3889*/
3890 if ( nind ) {
3891 t = pind;
3892 for ( i = 0; i < nind; i++ ) {
3893 r = t + 1;
3894 for ( j = i+1; j < nind; j++ ) {
3895 if ( *r < *t ) {
3896 k = *r; *r = *t; *t = k;
3897 }
3898 r++;
3899 }
3900 t++;
3901 }
3902 t = pind;
3903 *m++ = INDEX;
3904 *m++ = nind + 2;
3905 i = nind;
3906 NCOPY(m,t,i);
3907 }
3908/*
3909 #] Loose Vectors/Indices :
3910 #[ Vectors :
3911*/
3912 if ( nvec ) {
3913 t = pvec;
3914 for ( i = 2; i < nvec; i += 2 ) {
3915 r = t + 2;
3916 for ( j = i; j < nvec; j += 2 ) {
3917 if ( *r == *t ) {
3918 if ( *++r < t[1] ) {
3919 k = *r; *r = t[1]; t[1] = k;
3920 }
3921 r++;
3922 }
3923 else if ( *r < *t ) {
3924 k = *r; *r++ = *t; *t++ = k;
3925 k = *r; *r++ = *t; *t-- = k;
3926 }
3927 else { r += 2; }
3928 }
3929 t += 2;
3930 }
3931 t = pvec;
3932 *m++ = VECTOR;
3933 *m++ = nvec + 2;
3934 i = nvec;
3935 NCOPY(m,t,i);
3936 }
3937/*
3938 #] Vectors :
3939 #[ Dotproducts :
3940*/
3941 if ( ndot ) {
3942 to = m;
3943 m = t = pdot;
3944 i = ndot;
3945 while ( --i >= 0 ) {
3946 if ( *t > t[1] ) { j = *t; *t = t[1]; t[1] = j; }
3947 t += 3;
3948 }
3949 t = m;
3950 ndot *= 3;
3951 m += ndot;
3952 while ( t < (m-3) ) {
3953 r = t + 3;
3954 do {
3955 if ( *r == *t ) {
3956 if ( *++r == *++t ) {
3957 r++;
3958 if ( ( *r < MAXPOWER && t[1] < MAXPOWER )
3959 || ( *r > -MAXPOWER && t[1] > -MAXPOWER ) ) {
3960 t++;
3961 *t += *r;
3962 if ( *t > MAXPOWER || *t < -MAXPOWER ) {
3963 MLOCK(ErrorMessageLock);
3964 MesPrint("Exponent of dotproduct out of range: %d",*t);
3965 MUNLOCK(ErrorMessageLock);
3966 goto NormMin;
3967 }
3968 ndot -= 3;
3969 *r-- = *--m;
3970 *r-- = *--m;
3971 *r = *--m;
3972 if ( !*t ) {
3973 ndot -= 3;
3974 *t-- = *--m;
3975 *t-- = *--m;
3976 *t = *--m;
3977 t -= 3;
3978 break;
3979 }
3980 }
3981 else if ( *r < *++t ) {
3982 k = *r; *r++ = *t; *t = k;
3983 }
3984 else r++;
3985 t -= 2;
3986 }
3987 else if ( *r < *t ) {
3988 k = *r; *r++ = *t; *t++ = k;
3989 k = *r; *r++ = *t; *t = k;
3990 t -= 2;
3991 }
3992 else { r += 2; t--; }
3993 }
3994 else if ( *r < *t ) {
3995 k = *r; *r++ = *t; *t++ = k;
3996 k = *r; *r++ = *t; *t++ = k;
3997 k = *r; *r++ = *t; *t = k;
3998 t -= 2;
3999 }
4000 else { r += 3; }
4001 } while ( r < m );
4002 t += 3;
4003 }
4004 m = to;
4005 t = pdot;
4006 if ( ( i = ndot ) > 0 ) {
4007 *m++ = DOTPRODUCT;
4008 *m++ = i + 2;
4009 NCOPY(m,t,i);
4010 }
4011 }
4012/*
4013 #] Dotproducts :
4014 #[ Symbols :
4015*/
4016 if ( nsym ) {
4017 nsym <<= 1;
4018 t = psym;
4019 *m++ = SYMBOL;
4020 r = m;
4021 *m++ = ( i = nsym ) + 2;
4022 if ( i ) { do {
4023 if ( !*t ) {
4024 if ( t[1] < (2*MAXPOWER) ) { /* powers of i */
4025 if ( t[1] & 1 ) { *m++ = 0; *m++ = 1; }
4026 else *r -= 2;
4027 if ( *++t & 2 ) ncoef = -ncoef;
4028 t++;
4029 }
4030 }
4031 else if ( *t <= NumSymbols && *t > -2*MAXPOWER ) { /* Put powers in range */
4032 if ( ( ( ( t[1] > symbols[*t].maxpower ) && ( symbols[*t].maxpower < MAXPOWER ) ) ||
4033 ( ( t[1] < symbols[*t].minpower ) && ( symbols[*t].minpower > -MAXPOWER ) ) ) &&
4034 ( t[1] < 2*MAXPOWER ) && ( t[1] > -2*MAXPOWER ) ) {
4035 if ( i <= 2 || t[2] != *t ) goto NormZero;
4036 }
4037 if ( AN.ncmod == 1 && ( AC.modmode & ALSOPOWERS ) != 0 ) {
4038 if ( AC.cmod[0] == 1 ) t[1] = 0;
4039 else if ( t[1] >= 0 ) t[1] = 1 + (t[1]-1)%(AC.cmod[0]-1);
4040 else {
4041 t[1] = -1 - (-t[1]-1)%(AC.cmod[0]-1);
4042 if ( t[1] < 0 ) t[1] += (AC.cmod[0]-1);
4043 }
4044 }
4045 if ( ( t[1] < (2*MAXPOWER) && t[1] >= MAXPOWER )
4046 || ( t[1] > -(2*MAXPOWER) && t[1] <= -MAXPOWER ) ) {
4047 MLOCK(ErrorMessageLock);
4048 MesPrint("Exponent out of range: %d",t[1]);
4049 MUNLOCK(ErrorMessageLock);
4050 goto NormMin;
4051 }
4052 if ( AT.TrimPower && AR.PolyFunVar == *t && t[1] > AR.PolyFunPow ) {
4053 goto NormZero;
4054 }
4055 else if ( t[1] ) {
4056 *m++ = *t++;
4057 *m++ = *t++;
4058 }
4059 else { *r -= 2; t += 2; }
4060 }
4061 else {
4062 *m++ = *t++; *m++ = *t++;
4063 }
4064 } while ( (i-=2) > 0 ); }
4065 if ( *r <= 2 ) m = r-1;
4066 }
4067/*
4068 #] Symbols :
4069 #[ float_ :
4070
4071 Here we treat float_ functions and combined them with the regular
4072 coefficient.
4073*/
4074
4075#ifdef WITHFLOAT
4076 if ( withfloat ) {
4077 WORD floatsign = 3;
4078/*
4079 First check whether the coefficient is already 1/1
4080*/
4081 if ( ABS(ncoef) == 3 && n_coef[0] == 1 && n_coef[1] == 1 ) {
4082 if ( withfloat == 1 ) {
4083 t = firstfloat;
4084/*
4085 We only interfere by making the sign positive.
4086 The float_ function has already been tested.
4087*/
4088 if ( t[FUNHEAD+3] < 0 ) {
4089 t[FUNHEAD+3] = -t[FUNHEAD+3];
4090 floatsign = -floatsign;
4091 }
4092 if ( ncoef < 0 ) floatsign = -floatsign;
4093/*
4094 Copy to output;
4095*/
4096 AT.FloatPos = m-termout;
4097 i = t[1]; NCOPY(m,t,i)
4098 }
4099 else {
4100 PackFloat(m,aux4);
4101 if ( m[FUNHEAD+3] < 0 ) {
4102 m[FUNHEAD+3] = -m[FUNHEAD+3];
4103 floatsign = -floatsign;
4104 }
4105 if ( ncoef < 0 ) floatsign = -floatsign;
4106 AT.FloatPos = m-termout;
4107 m += m[1];
4108 }
4109 }
4110 else {
4111 if ( withfloat == 1 ) UnpackFloat(aux4,firstfloat);
4112 RatToFloat(aux5,(UWORD *)n_coef,ncoef);
4113 mpf_mul(aux4,aux4,aux5);
4114 PackFloat(m,aux4);
4115 AT.FloatPos = m-termout;
4116 if ( m[FUNHEAD+3] < 0 ) {
4117 m[FUNHEAD+3] = -m[FUNHEAD+3];
4118 floatsign = -floatsign;
4119 }
4120 m += m[1];
4121 }
4122 n_coef[0] = 1; n_coef[1] = 1; ncoef = floatsign;
4123 }
4124 else AT.FloatPos = 0;
4125#endif
4126
4127/*
4128 #] float_ :
4129 #[ Do Replace_ :
4130*/
4131 stop = (WORD *)(((UBYTE *)(termout)) + AM.MaxTer);
4132 i = ABS(ncoef);
4133 if ( ( m + i ) > stop ) {
4134 MLOCK(ErrorMessageLock);
4135 MesPrint("Term too complex during normalization");
4136 MUNLOCK(ErrorMessageLock);
4137 goto NormMin;
4138 }
4139 if ( AT.SS == AT.S0 ) {
4140 if ( ( m + i - termout ) > AT.SS->verbMaxTermSize ) {
4141 AT.SS->verbMaxTermSize = m+i-termout;
4142 }
4143 }
4144 if ( ReplaceType >= 0 ) {
4145 t = n_coef;
4146 i--;
4147 NCOPY(m,t,i);
4148 *m++ = ncoef;
4149 t = termout;
4150 *t = WORDDIF(m,t);
4151 if ( ReplaceType == 0 ) {
4152 AT.WorkPointer = termout+*termout;
4153 WildFill(BHEAD term,termout,AN.ReplaceScrat);
4154 termout = term + *term;
4155 }
4156 else {
4157 AT.WorkPointer = r = termout + *termout;
4158 WildFill(BHEAD r,termout,AN.ReplaceScrat);
4159 i = *r; m = term;
4160 NCOPY(m,r,i);
4161 termout = m;
4162
4163
4164 r = m = term;
4165 r += *term; r -= ABS(r[-1]);
4166 m++;
4167 while ( m < r ) {
4168 if ( *m >= FUNCTION && m[1] > FUNHEAD &&
4169 functions[*m-FUNCTION].spec != TENSORFUNCTION )
4170 m[2] |= DIRTYFLAG;
4171 m += m[1];
4172 }
4173 }
4174/*
4175 The next 'reset' cannot be done. We still need the expression
4176 in the buffer. Note though that this may cause a runaway pointer
4177 if we are not very careful.
4178
4179 C->numrhs = oldtoprhs;
4180 C->Pointer = C->Buffer + oldcpointer;
4181*/
4182 AT.NormDepth--;
4183 TermFree(n_llnum,"n_llnum");
4184 TermFree(n_coef,"NormCoef");
4185 return(1);
4186 }
4187 else {
4188 t = termout;
4189 k = WORDDIF(m,t);
4190 *t = k + i;
4191 m = term;
4192 NCOPY(m,t,k);
4193 i--;
4194 t = n_coef;
4195 NCOPY(m,t,i);
4196 *m++ = ncoef;
4197 }
4198/*
4199 #] Do Replace_ :
4200 #[ Errors and Finish :
4201*/
4202RegEnd:
4203 if ( termout < term + *term && termout >= term ) AT.WorkPointer = term + *term;
4204 else AT.WorkPointer = termout;
4205/*
4206 if ( termflag ) { We have to assign the term to $variable(s)
4207 TermAssign(term);
4208 }
4209*/
4210 AT.NormDepth--;
4211 TermFree(n_llnum,"n_llnum");
4212 TermFree(n_coef,"NormCoef");
4213 return(regval);
4214
4215NormInf:
4216 MLOCK(ErrorMessageLock);
4217 MesPrint("Division by zero during normalization");
4218 MUNLOCK(ErrorMessageLock);
4219 Terminate(-1);
4220
4221NormZZ:
4222 MLOCK(ErrorMessageLock);
4223 MesPrint("0^0 during normalization of term");
4224 MUNLOCK(ErrorMessageLock);
4225 Terminate(-1);
4226
4227NormPRF:
4228 MLOCK(ErrorMessageLock);
4229 MesPrint("0/0 in polyratfun during normalization of term");
4230 MUNLOCK(ErrorMessageLock);
4231 Terminate(-1);
4232
4233NormZero:
4234 *term = 0;
4235 AT.WorkPointer = termout;
4236 AT.NormDepth--;
4237 TermFree(n_llnum,"n_llnum");
4238 TermFree(n_coef,"NormCoef");
4239 return(regval);
4240
4241NormMin:
4242 AT.NormDepth--;
4243 TermFree(n_llnum,"n_llnum");
4244 TermFree(n_coef,"NormCoef");
4245 return(-1);
4246
4247FromNorm:
4248 MLOCK(ErrorMessageLock);
4249 MesCall("Norm");
4250 MUNLOCK(ErrorMessageLock);
4251 AT.NormDepth--;
4252 TermFree(n_llnum,"n_llnum");
4253 TermFree(n_coef,"NormCoef");
4254 return(-1);
4255
4256/*
4257 #] Errors and Finish :
4258*/
4259}
4260
4261/*
4262 #] Normalize :
4263 #[ ExtraSymbol :
4264*/
4265
4266int ExtraSymbol(WORD sym, WORD pow, WORD nsym, WORD *ppsym, WORD *ncoef)
4267{
4268 WORD *m, i;
4269 i = nsym;
4270 m = ppsym - 2;
4271 while ( i > 0 ) {
4272 if ( sym == *m ) {
4273 m++;
4274 if ( pow > 2*MAXPOWER || pow < -2*MAXPOWER
4275 || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
4276 MLOCK(ErrorMessageLock);
4277 MesPrint("Illegal wildcard power combination.");
4278 MUNLOCK(ErrorMessageLock);
4279 Terminate(-1);
4280 }
4281 *m += pow;
4282
4283 if ( ( sym <= NumSymbols && sym > -MAXPOWER )
4284 && ( symbols[sym].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
4285 *m %= symbols[sym].maxpower;
4286 if ( *m < 0 ) *m += symbols[sym].maxpower;
4287 if ( ( symbols[sym].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
4288 if ( ( ( symbols[sym].maxpower & 1 ) == 0 ) &&
4289 ( *m >= symbols[sym].maxpower/2 ) ) {
4290 *m -= symbols[sym].maxpower/2; *ncoef = -*ncoef;
4291 }
4292 }
4293 }
4294
4295 if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
4296 MLOCK(ErrorMessageLock);
4297 MesPrint("Power overflow during normalization");
4298 MUNLOCK(ErrorMessageLock);
4299 return(-1);
4300 }
4301 if ( !*m ) {
4302 m--;
4303 while ( i < nsym )
4304 { *m = m[2]; m++; *m = m[2]; m++; i++; }
4305 return(-1);
4306 }
4307 return(0);
4308 }
4309 else if ( sym < *m ) {
4310 m -= 2;
4311 i--;
4312 }
4313 else break;
4314 }
4315 m = ppsym;
4316 while ( i < nsym )
4317 { m--; m[2] = *m; m--; m[2] = *m; i++; }
4318 *m++ = sym;
4319 *m = pow;
4320 return(1);
4321}
4322
4323/*
4324 #] ExtraSymbol :
4325 #[ DoTheta :
4326*/
4327
4328int DoTheta(PHEAD WORD *t)
4329{
4330 GETBIDENTITY
4331 WORD k, *r1, *r2, *tstop, type;
4332 WORD ia, *ta, *tb, *stopa, *stopb;
4333 if ( AC.BracketNormalize ) return(-1);
4334 type = *t;
4335 k = t[1];
4336 tstop = t + k;
4337 t += FUNHEAD;
4338 if ( k <= FUNHEAD ) return(1);
4339 r1 = t;
4340 NEXTARG(r1)
4341 if ( r1 == tstop ) {
4342/*
4343 One argument
4344*/
4345 if ( *t == ARGHEAD ) {
4346 if ( type == THETA ) return(1);
4347 else return(0); /* THETA2 */
4348 }
4349 if ( *t < 0 ) {
4350 if ( *t == -SNUMBER ) {
4351 if ( t[1] < 0 ) return(0);
4352 else {
4353 if ( type == THETA2 && t[1] == 0 ) return(0);
4354 else return(1);
4355 }
4356 }
4357 return(-1);
4358 }
4359 k = t[*t-1];
4360 if ( *t == ABS(k)+1+ARGHEAD ) {
4361 if ( k > 0 ) return(1);
4362 else return(0);
4363 }
4364 return(-1);
4365 }
4366/*
4367 At least two arguments
4368*/
4369 r2 = r1;
4370 NEXTARG(r2)
4371 if ( r2 < tstop ) return(-1); /* More than 2 arguments */
4372/*
4373 Note now that zero has to be treated specially
4374 We take the criteria from the symmetrize routine
4375*/
4376 if ( *t == -SNUMBER && *r1 == -SNUMBER ) {
4377 if ( t[1] > r1[1] ) return(0);
4378 else if ( t[1] < r1[1] ) {
4379 return(1);
4380 }
4381 else if ( type == THETA ) return(1);
4382 else return(0); /* THETA2 */
4383 }
4384 else if ( t[1] == 0 && *t == -SNUMBER ) {
4385 if ( *r1 > 0 ) { }
4386 else if ( *t < *r1 ) return(1);
4387 else if ( *t > *r1 ) return(0);
4388 }
4389 else if ( r1[1] == 0 && *r1 == -SNUMBER ) {
4390 if ( *t > 0 ) { }
4391 else if ( *t < *r1 ) return(1);
4392 else if ( *t > *r1 ) return(0);
4393 }
4394 r2 = AT.WorkPointer;
4395 if ( *t < 0 ) {
4396 ta = r2;
4397 ToGeneral(t,ta,0);
4398 r2 += *r2;
4399 }
4400 else ta = t;
4401 if ( *r1 < 0 ) {
4402 tb = r2;
4403 ToGeneral(r1,tb,0);
4404 }
4405 else tb = r1;
4406 stopa = ta + *ta;
4407 stopb = tb + *tb;
4408 ta += ARGHEAD; tb += ARGHEAD;
4409 while ( ta < stopa ) {
4410 if ( tb >= stopb ) return(0);
4411 if ( ( ia = CompareTerms(BHEAD ta,tb,(WORD)1) ) < 0 ) return(0);
4412 if ( ia > 0 ) return(1);
4413 ta += *ta;
4414 tb += *tb;
4415 }
4416 if ( type == THETA ) return(1);
4417 else return(0); /* THETA2 */
4418}
4419
4420/*
4421 #] DoTheta :
4422 #[ DoDelta :
4423*/
4424
4425int DoDelta(WORD *t)
4426{
4427 WORD k, *r1, *r2, *tstop, isnum, isnum2, type = *t;
4428 if ( AC.BracketNormalize ) return(-1);
4429 k = t[1];
4430 if ( k <= FUNHEAD ) goto argzero;
4431 if ( k == FUNHEAD+ARGHEAD && t[FUNHEAD] == ARGHEAD ) goto argzero;
4432 tstop = t + k;
4433 t += FUNHEAD;
4434 r1 = t;
4435 NEXTARG(r1)
4436 if ( *t < 0 ) {
4437 k = 1;
4438 if ( *t == -SNUMBER ) { isnum = 1; k = t[1]; }
4439 else isnum = 0;
4440 }
4441 else {
4442 k = t[*t-1];
4443 k = ABS(k);
4444 if ( k == *t-ARGHEAD-1 ) isnum = 1;
4445 else isnum = 0;
4446 k = 1;
4447 }
4448 if ( r1 >= tstop ) { /* Single argument */
4449 if ( !isnum ) return(-1);
4450 if ( k == 0 ) goto argzero;
4451 goto argnonzero;
4452 }
4453 r2 = r1;
4454 NEXTARG(r2)
4455 if ( r2 < tstop ) return(-1);
4456 if ( *r1 < 0 ) {
4457 if ( *r1 == -SNUMBER ) { isnum2 = 1; }
4458 else isnum2 = 0;
4459 }
4460 else {
4461 k = r1[*r1-1];
4462 k = ABS(k);
4463 if ( k == *r1-ARGHEAD-1 ) isnum2 = 1;
4464 else isnum2 = 0;
4465 }
4466 if ( isnum != isnum2 ) return(-1);
4467 tstop = r1;
4468 while ( t < tstop && r1 < r2 ) {
4469 if ( *t != *r1 ) {
4470 if ( !isnum ) return(-1);
4471 goto argnonzero;
4472 }
4473 t++; r1++;
4474 }
4475 if ( t != tstop || r1 != r2 ) {
4476 if ( !isnum ) return(-1);
4477 goto argnonzero;
4478 }
4479argzero:
4480 if ( type == DELTA2 ) return(1);
4481 else return(0);
4482argnonzero:
4483 if ( type == DELTA2 ) return(0);
4484 else return(1);
4485}
4486
4487/*
4488 #] DoDelta :
4489 #[ DoRevert :
4490*/
4491
4492void DoRevert(WORD *fun, WORD *tmp)
4493{
4494 WORD *t, *r, *m, *to, *tt, *mm, i, j;
4495 to = fun + fun[1];
4496 r = fun + FUNHEAD;
4497 while ( r < to ) {
4498 if ( *r <= 0 ) {
4499 if ( *r == -REVERSEFUNCTION ) {
4500 m = r; mm = m+1;
4501 while ( mm < to ) *m++ = *mm++;
4502 to--;
4503 (fun[1])--;
4504 fun[2] |= DIRTYSYMFLAG;
4505 }
4506 else if ( *r <= -FUNCTION ) r++;
4507 else {
4508 if ( *r == -INDEX && r[1] < MINSPEC ) *r = -VECTOR;
4509 r += 2;
4510 }
4511 }
4512 else {
4513 if ( ( *r > ARGHEAD )
4514 && ( r[ARGHEAD+1] == REVERSEFUNCTION )
4515 && ( *r == (r[ARGHEAD]+ARGHEAD) )
4516 && ( r[ARGHEAD] == (r[ARGHEAD+2]+4) )
4517 && ( *(r+*r-3) == 1 )
4518 && ( *(r+*r-2) == 1 )
4519 && ( *(r+*r-1) == 3 ) ) {
4520 mm = r;
4521 r += ARGHEAD + 1;
4522 tt = r + r[1];
4523 r += FUNHEAD;
4524 m = tmp;
4525 t = r;
4526 j = 0;
4527 while ( t < tt ) {
4528 NEXTARG(t)
4529 j++;
4530 }
4531 while ( --j >= 0 ) {
4532 i = j;
4533 t = r;
4534 while ( --i >= 0 ) {
4535 NEXTARG(t)
4536 }
4537 if ( *t > 0 ) {
4538 i = *t;
4539 NCOPY(m,t,i);
4540 }
4541 else if ( *t <= -FUNCTION ) *m++ = *t++;
4542 else { *m++ = *t++; *m++ = *t++; }
4543 }
4544 i = WORDDIF(m,tmp);
4545 m = tmp;
4546 t = mm;
4547 r = t + *t;
4548 NCOPY(t,m,i);
4549 m = r;
4550 r = t;
4551 i = WORDDIF(to,m);
4552 NCOPY(t,m,i);
4553 fun[1] = WORDDIF(t,fun);
4554 to = t;
4555 fun[2] |= DIRTYSYMFLAG;
4556 }
4557 else r += *r;
4558 }
4559 }
4560}
4561
4562/*
4563 #] DoRevert :
4564 #] Normalize :
4565 #[ DetCommu :
4566
4567 Determines the number of terms in an expression that contain
4568 noncommuting objects. This can be used to see whether products of
4569 this expression can be evaluated with binomial coefficients.
4570
4571 We don't try to be fancy. If a term contains noncommuting objects
4572 we are not looking whether they can commute with complete other
4573 terms.
4574
4575 If the number gets too large we cut it off.
4576*/
4577
4578#define MAXNUMBEROFNONCOMTERMS 2
4579
4580WORD DetCommu(WORD *terms)
4581{
4582 WORD *t, *tnext, *tstop;
4583 WORD num = 0;
4584 if ( *terms == 0 ) return(0);
4585 if ( terms[*terms] == 0 ) return(0);
4586 t = terms;
4587 while ( *t ) {
4588 tnext = t + *t;
4589 tstop = tnext - ABS(tnext[-1]);
4590 t++;
4591 while ( t < tstop ) {
4592 if ( *t >= FUNCTION ) {
4593 if ( functions[*t-FUNCTION].commute ) {
4594 num++;
4595 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4596 break;
4597 }
4598 }
4599 else if ( *t == SUBEXPRESSION ) {
4600 if ( cbuf[t[4]].CanCommu[t[2]] ) {
4601 num++;
4602 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4603 break;
4604 }
4605 }
4606 else if ( *t == EXPRESSION ) {
4607 num++;
4608 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4609 break;
4610 }
4611 else if ( *t == DOLLAREXPRESSION ) {
4612/*
4613 Technically this is not correct. We have to test first
4614 whether this is DollarLocalCopy (in TFORM) and if so, use the
4615 local version. Anyway, this should be rare to never
4616 occurring because dollars should be replaced.
4617*/
4618 if ( cbuf[AM.dbufnum].CanCommu[t[2]] ) {
4619 num++;
4620 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4621 break;
4622 }
4623 }
4624 t += t[1];
4625 }
4626 t = tnext;
4627 }
4628 return(num);
4629}
4630
4631/*
4632 #] DetCommu :
4633 #[ DoesCommu :
4634
4635 Determines the number of noncommuting objects in a term.
4636 If the number gets too large we cut it off.
4637*/
4638
4639WORD DoesCommu(WORD *term)
4640{
4641 WORD *tstop;
4642 WORD num = 0;
4643 if ( *term == 0 ) return(0);
4644 tstop = term + *term;
4645 tstop = tstop - ABS(tstop[-1]);
4646 term++;
4647 while ( term < tstop ) {
4648 if ( ( *term >= FUNCTION ) && ( functions[*term-FUNCTION].commute ) ) {
4649 num++;
4650 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4651 }
4652 term += term[1];
4653 }
4654 return(num);
4655}
4656
4657/*
4658 #] DoesCommu :
4659 #[ PolyNormPoly :
4660
4661 Normalizes a polynomial
4662*/
4663
4664#ifdef EVALUATEGCD
4665WORD *PolyNormPoly (PHEAD WORD *Poly) {
4666
4667 GETBIDENTITY;
4668 WORD *buffer = AT.WorkPointer;
4669 WORD *p;
4670 if ( NewSort(BHEAD0) ) { Terminate(-1); }
4671 AR.CompareRoutine = (COMPAREDUMMY)(&CompareSymbols);
4672 while ( *Poly ) {
4673 p = Poly + *Poly;
4674 if ( SymbolNormalize(Poly) < 0 ) return(0);
4675 if ( StoreTerm(BHEAD Poly) ) {
4676 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
4677 LowerSortLevel();
4678 Terminate(-1);
4679 }
4680 Poly = p;
4681 }
4682 if ( EndSort(BHEAD buffer,1) < 0 ) {
4683 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
4684 Terminate(-1);
4685 }
4686 p = buffer;
4687 while ( *p ) p += *p;
4688 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
4689 AT.WorkPointer = p + 1;
4690 return(buffer);
4691}
4692#endif
4693
4694/*
4695 #] PolyNormPoly :
4696 #[ EvaluateGcd :
4697
4698 Try to evaluate the GCDFUNCTION gcd_.
4699 This function can have a number of arguments which can be numbers
4700 and/or polynomials. If there are objects that aren't SYMBOLS or numbers
4701 it cannot work currently.
4702
4703 To make this work properly we have to intervene in proces.c
4704 proces.c: if ( Normalize(BHEAD m) ) {
47051060 proces.c: if ( Normalize(BHEAD r) ) {
47061126?proces.c: if ( Normalize(BHEAD term) ) {
4707 proces.c: if ( Normalize(BHEAD AT.WorkPointer) ) goto PasErr;
47082308!proces.c: if ( ( retnorm = Normalize(BHEAD term) ) != 0 ) {
4709 proces.c: ReNumber(BHEAD term); Normalize(BHEAD term);
4710 proces.c: if ( Normalize(BHEAD v) ) Terminate(-1);
4711 proces.c: if ( Normalize(BHEAD w) ) { LowerSortLevel(); goto PolyCall; }
4712 proces.c: if ( Normalize(BHEAD term) ) goto PolyCall;
4713*/
4714#ifdef EVALUATEGCD
4715
4716WORD *EvaluateGcd(PHEAD WORD *subterm)
4717{
4718 GETBIDENTITY
4719 WORD *oldworkpointer = AT.WorkPointer, *work1, *work2, *work3;
4720 WORD *t, *tt, *ttt, *t1, *t2, *t3, *t4, *tstop;
4721 WORD ct, nnum;
4722 UWORD gcdnum, stor;
4723 WORD *lnum=n_llnum+1;
4724 WORD *num1, *num2, *num3, *den1, *den2, *den3;
4725 WORD sizenum1, sizenum2, sizenum3, sizeden1, sizeden2, sizeden3;
4726 int i, isnumeric = 0, numarg = 0 /*, sizearg */;
4727 LONG size;
4728/*
4729 Step 1: Look for -SNUMBER or -SYMBOL arguments.
4730 If encountered, treat everybody with it.
4731*/
4732 tt = subterm + subterm[1]; t = subterm + FUNHEAD;
4733
4734 while ( t < tt ) {
4735 numarg++;
4736 if ( *t == -SNUMBER ) {
4737 if ( t[1] == 0 ) {
4738gcdzero:;
4739 MLOCK(ErrorMessageLock);
4740 MesPrint("Trying to take the GCD involving a zero term.");
4741 MUNLOCK(ErrorMessageLock);
4742 return(0);
4743 }
4744 gcdnum = ABS(t[1]);
4745 t1 = subterm + FUNHEAD;
4746 while ( gcdnum > 1 && t1 < tt ) {
4747 if ( *t1 == -SNUMBER ) {
4748 stor = ABS(t1[1]);
4749 if ( stor == 0 ) goto gcdzero;
4750 if ( GcdLong(BHEAD (UWORD *)&stor,1,(UWORD *)&gcdnum,1,
4751 (UWORD *)lnum,&nnum) ) goto FromGCD;
4752 gcdnum = lnum[0];
4753 t1 += 2;
4754 continue;
4755 }
4756 else if ( *t1 == -SYMBOL ) goto gcdisone;
4757 else if ( *t1 < 0 ) goto gcdillegal;
4758/*
4759 Now we have to go through all the terms in the argument.
4760 This includes long numbers.
4761*/
4762 ttt = t1 + *t1;
4763 ct = *ttt; *ttt = 0;
4764 if ( t1[1] != 0 ) { /* First normalize the argument */
4765 t1 = PolyNormPoly(BHEAD t1+ARGHEAD);
4766 }
4767 else t1 += ARGHEAD;
4768 while ( *t1 ) {
4769 t1 += *t1;
4770 i = ABS(t1[-1]);
4771 t2 = t1 - i;
4772 i = (i-1)/2;
4773 t3 = t2+i-1;
4774 while ( t3 > t2 && *t3 == 0 ) { t3--; i--; }
4775 if ( GcdLong(BHEAD (UWORD *)t2,(WORD)i,(UWORD *)&gcdnum,1,
4776 (UWORD *)lnum,&nnum) ) {
4777 *ttt = ct;
4778 goto FromGCD;
4779 }
4780 gcdnum = lnum[0];
4781 if ( gcdnum == 1 ) {
4782 *ttt = ct;
4783 goto gcdisone;
4784 }
4785 }
4786 *ttt = ct;
4787 t1 = ttt;
4788 AT.WorkPointer = oldworkpointer;
4789 }
4790 if ( gcdnum == 1 ) goto gcdisone;
4791 oldworkpointer[0] = 4;
4792 oldworkpointer[1] = gcdnum;
4793 oldworkpointer[2] = 1;
4794 oldworkpointer[3] = 3;
4795 oldworkpointer[4] = 0;
4796 AT.WorkPointer = oldworkpointer + 5;
4797 return(oldworkpointer);
4798 }
4799 else if ( *t == -SYMBOL ) {
4800 t1 = subterm + FUNHEAD;
4801 i = t[1];
4802 while ( t1 < tt ) {
4803 if ( *t1 == -SNUMBER ) goto gcdisone;
4804 if ( *t1 == -SYMBOL ) {
4805 if ( t1[1] != i ) goto gcdisone;
4806 t1 += 2; continue;
4807 }
4808 if ( *t1 < 0 ) goto gcdillegal;
4809 ttt = t1 + *t1;
4810 ct = *ttt; *ttt = 0;
4811 if ( t1[1] != 0 ) { /* First normalize the argument */
4812 t2 = PolyNormPoly(BHEAD t1+ARGHEAD);
4813 }
4814 else t2 = t1 + ARGHEAD;
4815 while ( *t2 ) {
4816 t3 = t2+1;
4817 t2 = t2 + *t2;
4818 tstop = t2 - ABS(t2[-1]);
4819 while ( t3 < tstop ) {
4820 if ( *t3 != SYMBOL ) {
4821 *ttt = ct;
4822 goto gcdillegal;
4823 }
4824 t4 = t3 + 2;
4825 t3 += t3[1];
4826 while ( t4 < t3 ) {
4827 if ( *t4 == i && t4[1] > 0 ) goto nextterminarg;
4828 t4 += 2;
4829 }
4830 }
4831 *ttt = ct;
4832 goto gcdisone;
4833nextterminarg:;
4834 }
4835 *ttt = ct;
4836 t1 = ttt;
4837 AT.WorkPointer = oldworkpointer;
4838 }
4839 oldworkpointer[0] = 8;
4840 oldworkpointer[1] = SYMBOL;
4841 oldworkpointer[2] = 4;
4842 oldworkpointer[3] = t[1];
4843 oldworkpointer[4] = 1;
4844 oldworkpointer[5] = 1;
4845 oldworkpointer[6] = 1;
4846 oldworkpointer[7] = 3;
4847 oldworkpointer[8] = 0;
4848 AT.WorkPointer = oldworkpointer+9;
4849 return(oldworkpointer);
4850 }
4851 else if ( *t < 0 ) {
4852gcdillegal:;
4853 MLOCK(ErrorMessageLock);
4854 MesPrint("Illegal object in gcd_ function. Object not a number or a symbol.");
4855 MUNLOCK(ErrorMessageLock);
4856 goto FromGCD;
4857 }
4858 else if ( ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4859 else if ( t[1] != 0 ) {
4860 ttt = t + *t; ct = *ttt; *ttt = 0;
4861 t = PolyNormPoly(BHEAD t+ARGHEAD);
4862 *ttt = ct;
4863 if ( t[*t] == 0 && ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4864 AT.WorkPointer = oldworkpointer;
4865 t = ttt;
4866 }
4867 t += *t;
4868 }
4869/*
4870 At this point there are only generic arguments.
4871 There are however still two cases:
4872 1: There is an argument that is purely numerical
4873 In that case we have to take the gcd of all coefficients
4874 2: All arguments are nontrivial polynomials.
4875 Here we don't worry so much about the factor. (???)
4876 We know whether case 1 occurs when isnumeric > 0.
4877 We can look up numarg to get a good starting value.
4878*/
4879 AT.WorkPointer = oldworkpointer;
4880 if ( isnumeric ) {
4881 t = subterm + FUNHEAD;
4882 for ( i = 1; i < isnumeric; i++ ) {
4883 NEXTARG(t);
4884 }
4885 if ( t[1] != 0 ) { /* First normalize the argument */
4886 ttt = t + *t; ct = *ttt; *ttt = 0;
4887 t = PolyNormPoly(BHEAD t+ARGHEAD);
4888 *ttt = ct;
4889 }
4890 t += *t;
4891 i = (ABS(t[-1])-1)/2;
4892 den1 = t - 1 - i;
4893 num1 = den1 - i;
4894 sizenum1 = sizeden1 = i;
4895 while ( sizenum1 > 1 && num1[sizenum1-1] == 0 ) sizenum1--;
4896 while ( sizeden1 > 1 && den1[sizeden1-1] == 0 ) sizeden1--;
4897 work1 = AT.WorkPointer+1; work2 = work1+sizenum1;
4898 for ( i = 0; i < sizenum1; i++ ) work1[i] = num1[i];
4899 for ( i = 0; i < sizeden1; i++ ) work2[i] = den1[i];
4900 num1 = work1; den1 = work2;
4901 AT.WorkPointer = work2 = work2 + sizeden1;
4902 t = subterm + FUNHEAD;
4903 while ( t < tt ) {
4904 ttt = t + *t; ct = *ttt; *ttt = 0;
4905 if ( t[1] != 0 ) {
4906 t = PolyNormPoly(BHEAD t+ARGHEAD);
4907 }
4908 else t += ARGHEAD;
4909 while ( *t ) {
4910 t += *t;
4911 i = (ABS(t[-1])-1)/2;
4912 den2 = t - 1 - i;
4913 num2 = den2 - i;
4914 sizenum2 = sizeden2 = i;
4915 while ( sizenum2 > 1 && num2[sizenum2-1] == 0 ) sizenum2--;
4916 while ( sizeden2 > 1 && den2[sizeden2-1] == 0 ) sizeden2--;
4917 num3 = AT.WorkPointer;
4918 if ( GcdLong(BHEAD (UWORD *)num2,sizenum2,(UWORD *)num1,sizenum1,
4919 (UWORD *)num3,&sizenum3) ) goto FromGCD;
4920 sizenum1 = sizenum3;
4921 for ( i = 0; i < sizenum1; i++ ) num1[i] = num3[i];
4922 den3 = AT.WorkPointer;
4923 if ( GcdLong(BHEAD (UWORD *)den2,sizeden2,(UWORD *)den1,sizeden1,
4924 (UWORD *)den3,&sizeden3) ) goto FromGCD;
4925 sizeden1 = sizeden3;
4926 for ( i = 0; i < sizeden1; i++ ) den1[i] = den3[i];
4927 if ( sizenum1 == 1 && num1[0] == 1 && sizeden1 == 1 && den1[1] == 1 )
4928 goto gcdisone;
4929 }
4930 *ttt = ct;
4931 t = ttt;
4932 AT.WorkPointer = work2;
4933 }
4934 AT.WorkPointer = oldworkpointer;
4935/*
4936 Now copy the GCD to the 'output'
4937*/
4938 if ( sizenum1 > sizeden1 ) {
4939 while ( sizenum1 > sizeden1 ) den1[sizeden1++] = 0;
4940 }
4941 else if ( sizenum1 < sizeden1 ) {
4942 while ( sizenum1 < sizeden1 ) num1[sizenum1++] = 0;
4943 }
4944 t = oldworkpointer;
4945 i = 2*sizenum1+1;
4946 *t++ = i+1;
4947 if ( num1 != t ) { NCOPY(t,num1,sizenum1); }
4948 else t += sizenum1;
4949 if ( den1 != t ) { NCOPY(t,den1,sizeden1); }
4950 else t += sizeden1;
4951 *t++ = i;
4952 *t++ = 0;
4953 AT.WorkPointer = t;
4954 return(oldworkpointer);
4955 }
4956/*
4957 Now the real stuff with only polynomials.
4958 Pick up the shortest term to start.
4959 We are a bit brutish about this.
4960*/
4961 t = subterm + FUNHEAD;
4962 AT.WorkPointer += AM.MaxTer/sizeof(WORD);
4963 work2 = AT.WorkPointer;
4964/*
4965 sizearg = subterm[1];
4966*/
4967 i = 0; work3 = 0;
4968 while ( t < tt ) {
4969 i++;
4970 work1 = AT.WorkPointer;
4971 ttt = t + *t; ct = *ttt; *ttt = 0;
4972 t = PolyNormPoly(BHEAD t+ARGHEAD);
4973 if ( *work1 < AT.WorkPointer-work1 ) {
4974/*
4975 sizearg = AT.WorkPointer-work1;
4976*/
4977 numarg = i;
4978 work3 = work1;
4979 }
4980 *ttt = ct; t = ttt;
4981 }
4982 *AT.WorkPointer++ = 0;
4983/*
4984 We have properly normalized arguments and the shortest is indicated in work3
4985*/
4986 work1 = work3;
4987 while ( *work2 ) {
4988 if ( work2 != work3 ) {
4989 work1 = PolyGCD2(BHEAD work1,work2);
4990 }
4991 while ( *work2 ) work2 += *work2;
4992 work2++;
4993 }
4994 work2 = work1;
4995 while ( *work2 ) work2 += *work2;
4996 size = work2 - work1 + 1;
4997 t = oldworkpointer;
4998 NCOPY(t,work1,size);
4999 AT.WorkPointer = t;
5000 return(oldworkpointer);
5001
5002gcdisone:;
5003 oldworkpointer[0] = 4;
5004 oldworkpointer[1] = 1;
5005 oldworkpointer[2] = 1;
5006 oldworkpointer[3] = 3;
5007 oldworkpointer[4] = 0;
5008 AT.WorkPointer = oldworkpointer+5;
5009 return(oldworkpointer);
5010FromGCD:
5011 MLOCK(ErrorMessageLock);
5012 MesCall("EvaluateGcd");
5013 MUNLOCK(ErrorMessageLock);
5014 return(0);
5015}
5016
5017#endif
5018
5019/*
5020 #] EvaluateGcd :
5021 #[ TreatPolyRatFun :
5022
5023 if ( AR.PolyFunExp == 1 ) we have to trim the contents of the polyratfun
5024 down to its most divergent term and give it coefficient +1. This is done
5025 by taking the terms with the least power in the variable in the numerator
5026 and in the denominator and then combine them.
5027 Answer is either PolyRatFun(ep^n,1) or PolyRatFun(1,1) or PolyRatFun(1,ep^n)
5028*/
5029
5030int TreatPolyRatFun(PHEAD WORD *prf)
5031{
5032 WORD *t, *tstop, *r, *rstop, *m, *mstop;
5033 WORD exp1 = MAXPOWER, exp2 = MAXPOWER;
5034 t = prf+FUNHEAD;
5035 if ( *t < 0 ) {
5036 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
5037 if ( exp1 > 1 ) exp1 = 1;
5038 t += 2;
5039 }
5040 else {
5041 if ( exp1 > 0 ) exp1 = 0;
5042 NEXTARG(t)
5043 }
5044 }
5045 else {
5046 tstop = t + *t;
5047 t += ARGHEAD;
5048 while ( t < tstop ) {
5049/*
5050 Now look for the minimum power of AR.PolyFunVar
5051*/
5052 r = t+1;
5053 t += *t;
5054 rstop = t - ABS(t[-1]);
5055 while ( r < rstop ) {
5056 if ( *r != SYMBOL ) { r += r[1]; continue; }
5057 m = r;
5058 mstop = m + m[1];
5059 m += 2;
5060 while ( m < mstop ) {
5061 if ( *m == AR.PolyFunVar ) {
5062 if ( m[1] < exp1 ) exp1 = m[1];
5063 break;
5064 }
5065 m += 2;
5066 }
5067 if ( m == mstop ) {
5068 if ( exp1 > 0 ) exp1 = 0;
5069 }
5070 break;
5071 }
5072 if ( r == rstop ) {
5073 if ( exp1 > 0 ) exp1 = 0;
5074 }
5075 }
5076 t = tstop;
5077 }
5078 if ( *t < 0 ) {
5079 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
5080 if ( exp2 > 1 ) exp2 = 1;
5081 }
5082 else {
5083 if ( exp2 > 0 ) exp2 = 0;
5084 }
5085 }
5086 else {
5087 tstop = t + *t;
5088 t += ARGHEAD;
5089 while ( t < tstop ) {
5090/*
5091 Now look for the minimum power of AR.PolyFunVar
5092*/
5093 r = t+1;
5094 t += *t;
5095 rstop = t - ABS(t[-1]);
5096 while ( r < rstop ) {
5097 if ( *r != SYMBOL ) { r += r[1]; continue; }
5098 m = r;
5099 mstop = m + m[1];
5100 m += 2;
5101 while ( m < mstop ) {
5102 if ( *m == AR.PolyFunVar ) {
5103 if ( m[1] < exp2 ) exp2 = m[1];
5104 break;
5105 }
5106 m += 2;
5107 }
5108 if ( m == mstop ) {
5109 if ( exp2 > 0 ) exp2 = 0;
5110 }
5111 break;
5112 }
5113 if ( r == rstop ) {
5114 if ( exp2 > 0 ) exp2 = 0;
5115 }
5116 }
5117 }
5118/*
5119 Now we can compose the output.
5120 Notice that the output can never be longer than the input provided
5121 we never can have arguments that consist of just a function.
5122*/
5123 exp1 = exp1-exp2;
5124/* if ( exp1 > 0 ) exp1 = 0; */
5125 t = prf+FUNHEAD;
5126 if ( exp1 == 0 ) {
5127 *t++ = -SNUMBER; *t++ = 1;
5128 *t++ = -SNUMBER; *t++ = 1;
5129 }
5130 else if ( exp1 > 0 ) {
5131 if ( exp1 == 1 ) {
5132 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
5133 }
5134 else {
5135 *t++ = 8+ARGHEAD;
5136 *t++ = 0;
5137 FILLARG(t);
5138 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
5139 *t++ = exp1; *t++ = 1; *t++ = 1; *t++ = 3;
5140 }
5141 *t++ = -SNUMBER; *t++ = 1;
5142 }
5143 else {
5144 *t++ = -SNUMBER; *t++ = 1;
5145 if ( exp1 == -1 ) {
5146 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
5147 }
5148 else {
5149 *t++ = 8+ARGHEAD;
5150 *t++ = 0;
5151 FILLARG(t);
5152 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
5153 *t++ = -exp1; *t++ = 1; *t++ = 1; *t++ = 3;
5154 }
5155 }
5156 prf[2] = 0; /* Clean */
5157 prf[1] = t - prf;
5158 return(0);
5159}
5160
5161/*
5162 #] TreatPolyRatFun :
5163 #[ DropCoefficient :
5164*/
5165
5166void DropCoefficient(PHEAD WORD *term)
5167{
5168 GETBIDENTITY
5169 WORD *t = term + *term;
5170 WORD n, na;
5171 n = t[-1]; na = ABS(n);
5172 t -= na;
5173 if ( n == 3 && t[0] == 1 && t[1] == 1 ) return;
5174 *AN.RepPoint = 1;
5175 t[0] = 1; t[1] = 1; t[2] = 3;
5176 *term -= (na-3);
5177}
5178
5179/*
5180 #] DropCoefficient :
5181 #[ DropSymbols :
5182*/
5183
5184void DropSymbols(PHEAD WORD *term)
5185{
5186 GETBIDENTITY
5187 WORD *tend = term + *term, *t1, *t2, *tstop;
5188 tstop = tend - ABS(tend[-1]);
5189 t1 = term+1;
5190 while ( t1 < tstop ) {
5191 if ( *t1 == SYMBOL ) {
5192 *AN.RepPoint = 1;
5193 t2 = t1+t1[1];
5194 while ( t2 < tend ) *t1++ = *t2++;
5195 *term = t1 - term;
5196 break;
5197 }
5198 t1 += t1[1];
5199 }
5200}
5201
5202/*
5203 #] DropSymbols :
5204 #[ SymbolNormalize :
5205*/
5214int SymbolNormalize(WORD *term)
5215{
5216 GETIDENTITY
5217 WORD *t, *b, *bb, *tt, *m, *tstop;
5218 // Here we use a stack-allocated array, since things are much smaller
5219 // compared to the full Normalize routine.
5220 WORD buffer[7*NORMSIZE];
5221 int i;
5222 b = buffer;
5223 *b++ = SYMBOL; *b++ = 2;
5224 t = term + *term;
5225 tstop = t - ABS(t[-1]);
5226 t = term + 1;
5227 while ( t < tstop ) { /* Step 1: collect symbols */
5228 if ( *t == SYMBOL && t < tstop ) {
5229 for ( i = 2; i < t[1]; i += 2 ) {
5230 const WORD sym = t[i];
5231 const WORD pow = t[i+1];
5232 bb = buffer+2;
5233 while ( bb < b ) {
5234 if ( bb[0] == sym ) { /* add powers */
5235 bb[1] += pow;
5236 if ( bb[1] > MAXPOWER || bb[1] < -MAXPOWER ) {
5237 MLOCK(ErrorMessageLock);
5238 MesPrint("Power in SymbolNormalize out of range");
5239 MUNLOCK(ErrorMessageLock);
5240 return(-1);
5241 }
5242 if ( bb[1] == 0 ) {
5243 b -= 2;
5244 while ( bb < b ) {
5245 bb[0] = bb[2]; bb[1] = bb[3]; bb += 2;
5246 }
5247 }
5248 goto Nexti;
5249 }
5250 else if ( bb[0] > sym ) { /* insert it */
5251 m = b;
5252 while ( m > bb ) { m[1] = m[-1]; m[0] = m[-2]; m -= 2; }
5253 b += 2;
5254 bb[0] = sym;
5255 bb[1] = pow;
5256 goto Nexti;
5257 }
5258 bb += 2;
5259 }
5260 if ( bb >= b ) { /* add it to the end */
5261 *b++ = sym; *b++ = pow;
5262 }
5263Nexti:;
5264 }
5265 }
5266 else {
5267 MLOCK(ErrorMessageLock);
5268 MesPrint("Illegal term in SymbolNormalize");
5269 MUNLOCK(ErrorMessageLock);
5270 return(-1);
5271 }
5272 t += t[1];
5273 }
5274 buffer[1] = b - buffer;
5275/*
5276 Veto negative powers
5277*/
5278 if ( AT.LeaveNegative == 0 ) {
5279 b = buffer; bb = b + b[1]; b += 3;
5280 while ( b < bb ) {
5281 if ( *b < 0 ) {
5282 MLOCK(ErrorMessageLock);
5283 MesPrint("Negative power in SymbolNormalize");
5284 MUNLOCK(ErrorMessageLock);
5285 return(-1);
5286 }
5287 b += 2;
5288 }
5289 }
5290/*
5291 Now we use the fact that the new term will not be longer than the old one
5292 Actually it should be shorter when there is more than one subterm!
5293 Copy back.
5294*/
5295 i = buffer[1];
5296 b = buffer; tt = term + 1;
5297 if ( i > 2 ) { NCOPY(tt,b,i) }
5298 if ( tt < tstop ) {
5299 i = term[*term-1];
5300 if ( i < 0 ) i = -i;
5301 *term -= (tstop-tt);
5302 NCOPY(tt,tstop,i)
5303 }
5304 return(0);
5305}
5306
5307/*
5308 #] SymbolNormalize :
5309 #[ TestFunFlag :
5310
5311 Tests whether a function still has unsubstituted subexpressions
5312 This function has its dirtyflag on!
5313*/
5314
5315int TestFunFlag(PHEAD WORD *tfun)
5316{
5317 WORD *t, *tstop, *r, *rstop, *m, *mstop;
5318 /* DENOMINATOR and EXPONENT functions have spec <= 0, but may have
5319 * subexpression arguments to search for. */
5320 if ( *tfun != DENOMINATOR && *tfun != EXPONENT
5321 && functions[*tfun-FUNCTION].spec <= 0 ) {
5322
5323 return(0);
5324 }
5325 tstop = tfun + tfun[1];
5326 t = tfun + FUNHEAD;
5327 while ( t < tstop ) {
5328 if ( *t < 0 ) { NEXTARG(t); continue; }
5329 rstop = t + *t;
5330 if ( t[1] == 0 ) { t = rstop; continue; }
5331 r = t + ARGHEAD;
5332 while ( r < rstop ) { /* Here we loop over terms */
5333 m = r+1; mstop = r+*r; mstop -= ABS(mstop[-1]);
5334 while ( m < mstop ) { /* Loop over the subterms */
5335 if ( *m == SUBEXPRESSION || *m == EXPRESSION || *m == DOLLAREXPRESSION ) return(1);
5336 if ( ( *m >= FUNCTION ) && ( ( m[2] & DIRTYFLAG ) == DIRTYFLAG )
5337 && ( *m != REPLACEMENT ) && TestFunFlag(BHEAD m) ) return(1);
5338 m += m[1];
5339 }
5340 r += *r;
5341 }
5342 t += *t;
5343 }
5344 return(0);
5345}
5346
5347/*
5348 #] TestFunFlag :
5349 #[ BracketNormalize :
5350*/
5351
5352int BracketNormalize(PHEAD WORD *term)
5353{
5354 WORD *stop = term+*term-3, *t, *tt, *tstart, *r;
5355 WORD *oldwork = AT.WorkPointer;
5356 WORD *termout;
5357 WORD i, ii, j;
5358 termout = AT.WorkPointer = term+*term;
5359/*
5360 First collect all functions and sort them
5361*/
5362 tt = termout+1; t = term+1;
5363 while ( t < stop ) {
5364 if ( *t >= FUNCTION ) { i = t[1]; NCOPY(tt,t,i); }
5365 else t += t[1];
5366 }
5367 if ( tt > termout+1 && tt-termout-1 > termout[2] ) { /* sorting */
5368 r = termout+1; ii = tt-r;
5369 for ( i = 0; i < ii-FUNHEAD; i += FUNHEAD ) { /* Bubble sort */
5370 for ( j = i+FUNHEAD; j > 0; j -= FUNHEAD ) {
5371 if ( functions[r[j-FUNHEAD]-FUNCTION].commute
5372 && functions[r[j]-FUNCTION].commute == 0 ) break;
5373 if ( r[j-FUNHEAD] > r[j] ) EXCH(r[j-FUNHEAD],r[j])
5374 else break;
5375 }
5376 }
5377 }
5378
5379 tstart = tt; t = term + 1; *tt++ = DELTA; *tt++ = 2;
5380 while ( t < stop ) {
5381 if ( *t == DELTA ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5382 else t += t[1];
5383 }
5384 if ( tstart[1] > 2 ) {
5385 for ( r = tstart+2; r < tstart+tstart[1]; r += 2 ) {
5386 if ( r[0] > r[1] ) EXCH(r[0],r[1])
5387 }
5388 }
5389 if ( tstart[1] > 4 ) { /* sorting */
5390 r = tstart+2; ii = tstart[1]-2;
5391 for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5392 for ( j = i+2; j > 0; j -= 2 ) {
5393 if ( r[j-2] > r[j] ) {
5394 EXCH(r[j-2],r[j])
5395 EXCH(r[j-1],r[j+1])
5396 }
5397 else if ( r[j-2] < r[j] ) break;
5398 else {
5399 if ( r[j-1] > r[j+1] ) EXCH(r[j-1],r[j+1])
5400 else break;
5401 }
5402 }
5403 }
5404 tt = tstart+tstart[1];
5405 }
5406 else if ( tstart[1] == 2 ) { tt = tstart; }
5407 else tt = tstart+4;
5408
5409 tstart = tt; t = term + 1; *tt++ = INDEX; *tt++ = 2;
5410 while ( t < stop ) {
5411 if ( *t == INDEX ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5412 else t += t[1];
5413 }
5414 if ( tstart[1] >= 4 ) { /* sorting */
5415 r = tstart+2; ii = tstart[1]-2;
5416 for ( i = 0; i < ii-1; i += 1 ) { /* Bubble sort */
5417 for ( j = i+1; j > 0; j -= 1 ) {
5418 if ( r[j-1] > r[j] ) EXCH(r[j-1],r[j])
5419 else break;
5420 }
5421 }
5422 tt = tstart+tstart[1];
5423 }
5424 else if ( tstart[1] == 2 ) { tt = tstart; }
5425 else tt = tstart+3;
5426
5427 tstart = tt; t = term + 1; *tt++ = DOTPRODUCT; *tt++ = 2;
5428 while ( t < stop ) {
5429 if ( *t == DOTPRODUCT ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5430 else t += t[1];
5431 }
5432 if ( tstart[1] > 5 ) { /* sorting */
5433 r = tstart+2; ii = tstart[1]-2;
5434 for ( i = 0; i < ii; i += 3 ) {
5435 if ( r[i] > r[i+1] ) EXCH(r[i],r[i+1])
5436 }
5437 for ( i = 0; i < ii-3; i += 3 ) { /* Bubble sort */
5438 for ( j = i+3; j > 0; j -= 3 ) {
5439 if ( r[j-3] < r[j] ) break;
5440 if ( r[j-3] > r[j] ) {
5441 EXCH(r[j-3],r[j])
5442 EXCH(r[j-2],r[j+1])
5443 }
5444 else {
5445 if ( r[j-2] > r[j+1] ) EXCH(r[j-2],r[j+1])
5446 else break;
5447 }
5448 }
5449 }
5450 tt = tstart+tstart[1];
5451 }
5452 else if ( tstart[1] == 2 ) { tt = tstart; }
5453 else {
5454 if ( tstart[2] > tstart[3] ) EXCH(tstart[2],tstart[3])
5455 tt = tstart+5;
5456 }
5457
5458 tstart = tt; t = term + 1; *tt++ = SYMBOL; *tt++ = 2;
5459 while ( t < stop ) {
5460 if ( *t == SYMBOL ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5461 else t += t[1];
5462 }
5463 if ( tstart[1] > 4 ) { /* sorting */
5464 r = tstart+2; ii = tstart[1]-2;
5465 for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5466 for ( j = i+2; j > 0; j -= 2 ) {
5467 if ( r[j-2] > r[j] ) EXCH(r[j-2],r[j])
5468 else break;
5469 }
5470 }
5471 tt = tstart+tstart[1];
5472 }
5473 else if ( tstart[1] == 2 ) { tt = tstart; }
5474 else tt = tstart+4;
5475
5476 tstart = tt; t = term + 1; *tt++ = SETSET; *tt++ = 2;
5477 while ( t < stop ) {
5478 if ( *t == SETSET ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5479 else t += t[1];
5480 }
5481 if ( tstart[1] > 4 ) { /* sorting */
5482 r = tstart+2; ii = tstart[1]-2;
5483 for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5484 for ( j = i+2; j > 0; j -= 2 ) {
5485 if ( r[j-2] > r[j] ) {
5486 EXCH(r[j-2],r[j])
5487 EXCH(r[j-1],r[j+1])
5488 }
5489 else break;
5490 }
5491 }
5492 tt = tstart+tstart[1];
5493 }
5494 else if ( tstart[1] == 2 ) { tt = tstart; }
5495 else tt = tstart+4;
5496 *tt++ = 1; *tt++ = 1; *tt++ = 3;
5497 t = term; i = *termout = tt - termout; tt = termout;
5498 NCOPY(t,tt,i);
5499 AT.WorkPointer = oldwork;
5500 return(0);
5501}
5502
5503/*
5504 #] BracketNormalize :
5505*/
AddRHS
WORD * AddRHS(int num, int type)
Definition comtool.c:210
DoubleCbuffer
WORD * DoubleCbuffer(int num, WORD *w, int par)
Definition comtool.c:143
GetFirstTerm
int GetFirstTerm(WORD *, int, int)
Definition execute.c:1975
CompCoef
WORD CompCoef(WORD *, WORD *)
Definition reken.c:3066
EndSort
LONG EndSort(PHEAD WORD *, int)
Definition sort.c:488
LowerSortLevel
void LowerSortLevel(void)
Definition sort.c:4731
StoreTerm
int StoreTerm(PHEAD WORD *)
Definition sort.c:4311
NextPrime
WORD NextPrime(PHEAD WORD)
Definition reken.c:3685
NewSort
int NewSort(PHEAD0)
Definition sort.c:397
Compare1
WORD Compare1(PHEAD WORD *, WORD *, WORD)
Definition sort.c:2393
CompareSymbols
WORD CompareSymbols(PHEAD WORD *, WORD *, WORD)
Definition sort.c:2856
SymbolNormalize
int SymbolNormalize(WORD *term)
Definition normal.c:5214
CbUf::Top
WORD * Top
Definition structs.h:972
CbUf::rhs
WORD ** rhs
Definition structs.h:975
CbUf::Buffer
WORD * Buffer
Definition structs.h:971
CbUf::Pointer
WORD * Pointer
Definition structs.h:973
Generated by doxygen 1.9.8