FORM v5.0.1-23-g7a8f756
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 MLOCK(ErrorMessageLock);
2512 MesPrint("Illegal code in Norm");
2513#ifdef DEBUGON
2514 {
2515 UBYTE OutBuf[140];
2516 AO.OutFill = AO.OutputLine = OutBuf;
2517 t = term;
2518 AO.OutSkip = 3;
2519 FiniLine();
2520 i = *t;
2521 while ( --i >= 0 ) {
2522 TalToLine((UWORD)(*t++));
2523 TokenToLine((UBYTE *)" ");
2524 }
2525 AO.OutSkip = 0;
2526 FiniLine();
2527 }
2528#endif
2529 MUNLOCK(ErrorMessageLock);
2530 goto NormMin;
2531 }
2532 if ( *t == REPLACEMENT ) {
2533 if ( AR.Eside != LHSIDE ) ReplaceVeto--;
2534 pcom[ncom++] = t;
2535 break;
2536 }
2537/*
2538 if ( *t == AM.termfunnum && t[1] == FUNHEAD+2
2539 && t[FUNHEAD] == -DOLLAREXPRESSION ) termflag++;
2540*/
2541 if ( *t == DUMMYFUN || *t == DUMMYTEN ) {}
2542 else {
2543 if ( *t < (FUNCTION + WILDOFFSET) ) {
2544 if ( ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2545 || ( functions[*t-FUNCTION].minnumargs > 0 ) )
2546 && ( ( t[2] & DIRTYFLAG ) != 0 ) ) {
2547/*
2548 Number of arguments is bounded. And we have not checked.
2549*/
2550 WORD *ta = t + FUNHEAD, *tb = t + t[1];
2551 int numarg = 0;
2552 while ( ta < tb ) { numarg++; NEXTARG(ta) }
2553 if ( ( functions[*t-FUNCTION].maxnumargs > 0 )
2554 && ( numarg >= functions[*t-FUNCTION].maxnumargs ) )
2555 goto NormZero;
2556 if ( ( functions[*t-FUNCTION].minnumargs > 0 )
2557 && ( numarg < functions[*t-FUNCTION].minnumargs ) )
2558 goto NormZero;
2559 }
2560doflags:
2561 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION].tabl == 0 ) ) {
2562 t[2] &= ~DIRTYFLAG;
2563 t[2] |= DIRTYSYMFLAG;
2564 }
2565 if ( functions[*t-FUNCTION].commute ) { pnco[nnco++] = t; }
2566 else { pcom[ncom++] = t; }
2567 }
2568 else {
2569 if ( ( ( t[2] & DIRTYFLAG ) != 0 ) && ( functions[*t-FUNCTION-WILDOFFSET].tabl == 0 ) ) {
2570 t[2] &= ~DIRTYFLAG;
2571 t[2] |= DIRTYSYMFLAG;
2572 }
2573 if ( functions[*t-FUNCTION-WILDOFFSET].commute ) {
2574 pnco[nnco++] = t;
2575 }
2576 else { pcom[ncom++] = t; }
2577 }
2578 }
2579
2580 /* Now hunt for contractible indices */
2581
2582 if ( ( *t < (FUNCTION + WILDOFFSET)
2583 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) || (
2584 *t >= (FUNCTION + WILDOFFSET)
2585 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION ) ) {
2586 if ( *t >= GAMMA && *t <= GAMMASEVEN ) t++;
2587 t += FUNHEAD;
2588 while ( t < r ) {
2589 if ( *t == AM.vectorzero ) goto NormZero;
2590 if ( *t >= AM.OffsetIndex && ( *t >= AM.DumInd
2591 || ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2592 pcon[ncon++] = t;
2593 }
2594 else if ( *t == FUNNYWILD ) { t++; }
2595 t++;
2596 }
2597 }
2598 else {
2599 t += FUNHEAD;
2600 while ( t < r ) {
2601 if ( *t > 0 ) {
2602/*
2603 Here we should worry about a recursion
2604 A problem is the possibility of a construct
2605 like f(mu+nu)
2606*/
2607 t += *t;
2608 }
2609 else if ( *t <= -FUNCTION ) t++;
2610 else if ( *t == -INDEX ) {
2611 if ( t[1] >= AM.OffsetIndex &&
2612 ( t[1] >= AM.DumInd || ( t[1] < AM.WilInd
2613 && indices[t[1]-AM.OffsetIndex].dimension ) ) )
2614 pcon[ncon++] = t+1;
2615 t += 2;
2616 }
2617 else if ( *t == -SYMBOL ) {
2618 if ( t[1] >= MAXPOWER && t[1] < 2*MAXPOWER ) {
2619 *t = -SNUMBER;
2620 t[1] -= MAXPOWER;
2621 }
2622 else if ( t[1] < -MAXPOWER && t[1] > -2*MAXPOWER ) {
2623 *t = -SNUMBER;
2624 t[1] += MAXPOWER;
2625 }
2626 else t += 2;
2627 }
2628 else t += 2;
2629 }
2630 }
2631 break;
2632 }
2633 t = r;
2634TryAgain:;
2635 } while ( t < m );
2636 if ( ANsc ) {
2637 AN.cTerm = ANsc;
2638 r = t = ANsr; m = ANsm;
2639 ANsc = ANsm = ANsr = 0;
2640 goto conscan;
2641 }
2642/*
2643 #] First scan :
2644 #[ Easy denominators :
2645
2646 Easy denominators are denominators that can be replaced by
2647 negative powers of individual subterms. This may add to all
2648 our sublists.
2649
2650*/
2651 if ( nden ) {
2652 for ( k = 0, i = 0; i < nden; i++ ) {
2653 t = pden[i];
2654 if ( ( t[2] & DIRTYFLAG ) == 0 ) continue;
2655 r = t + t[1]; m = t + FUNHEAD;
2656 if ( m >= r ) {
2657 for ( j = i+1; j < nden; j++ ) pden[j-1] = pden[j];
2658 nden--;
2659 for ( j = 0; j < nnco; j++ ) if ( pnco[j] == t ) break;
2660 for ( j++; j < nnco; j++ ) pnco[j-1] = pnco[j];
2661 nnco--;
2662 i--;
2663 }
2664 else {
2665 NEXTARG(m);
2666 if ( m >= r ) continue;
2667/*
2668 We have more than one argument. Split the function.
2669*/
2670 if ( k == 0 ) {
2671 k = 1; to = termout; from = term;
2672 }
2673 while ( from < t ) *to++ = *from++;
2674 m = t + FUNHEAD;
2675 while ( m < r ) {
2676 stop = to;
2677 *to++ = DENOMINATOR;
2678 for ( j = 1; j < FUNHEAD; j++ ) *to++ = 0;
2679 if ( *m < -FUNCTION ) *to++ = *m++;
2680 else if ( *m < 0 ) { *to++ = *m++; *to++ = *m++; }
2681 else {
2682 j = *m; while ( --j >= 0 ) *to++ = *m++;
2683 }
2684 stop[1] = WORDDIF(to,stop);
2685 }
2686 from = r;
2687 if ( i == nden - 1 ) {
2688 stop = term + *term;
2689 while ( from < stop ) *to++ = *from++;
2690 i = *termout = WORDDIF(to,termout);
2691 to = term; from = termout;
2692 while ( --i >= 0 ) *to++ = *from++;
2693 goto Restart;
2694 }
2695 }
2696 }
2697 for ( i = 0; i < nden; i++ ) {
2698 t = pden[i];
2699 if ( ( t[2] & DIRTYFLAG ) == 0 ) continue;
2700 t[2] = 0;
2701 if ( t[FUNHEAD] == -SYMBOL ) {
2702 WORD change;
2703 t += FUNHEAD+1;
2704 change = ExtraSymbol(*t,-1,nsym,ppsym,&ncoef);
2705 nsym += change;
2706 ppsym += change * 2;
2707 goto DropDen;
2708 }
2709 else if ( t[FUNHEAD] == -SNUMBER ) {
2710 t += FUNHEAD+1;
2711 if ( *t == 0 ) goto NormInf;
2712 if ( *t < 0 ) { *AT.WorkPointer = -*t; j = -1; }
2713 else { *AT.WorkPointer = *t; j = 1; }
2714 ncoef = REDLENG(ncoef);
2715 if ( Divvy(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)AT.WorkPointer,j) )
2716 goto FromNorm;
2717 ncoef = INCLENG(ncoef);
2718 goto DropDen;
2719 }
2720 else if ( t[FUNHEAD] == ARGHEAD ) goto NormInf;
2721 else if ( t[FUNHEAD] > 0 && t[FUNHEAD+ARGHEAD] ==
2722 t[FUNHEAD]-ARGHEAD ) {
2723 /* Only one term */
2724 r = t + t[1] - 1;
2725 t += FUNHEAD + ARGHEAD + 1;
2726 j = *r;
2727 m = r - ABS(*r) + 1;
2728 if ( j != 3 || ( ( *m != 1 ) || ( m[1] != 1 ) ) ) {
2729 ncoef = REDLENG(ncoef);
2730 if ( DivRat(BHEAD (UWORD *)n_coef,ncoef,(UWORD *)m,REDLENG(j),(UWORD *)n_coef,&ncoef) ) goto FromNorm;
2731 ncoef = INCLENG(ncoef);
2732 j = ABS(j) - 3;
2733 t[-FUNHEAD-ARGHEAD] -= j;
2734 t[-ARGHEAD-1] -= j;
2735 t[-1] -= j;
2736 m[0] = m[1] = 1;
2737 m[2] = 3;
2738 }
2739 while ( t < m ) {
2740 r = t + t[1];
2741 if ( *t == SYMBOL || *t == DOTPRODUCT ) {
2742 k = t[1];
2743 pden[i][1] -= k;
2744 pden[i][FUNHEAD] -= k;
2745 pden[i][FUNHEAD+ARGHEAD] -= k;
2746 m -= k;
2747 stop = m + 3;
2748 tt = to = t;
2749 from = r;
2750 if ( *t == SYMBOL ) {
2751 t += 2;
2752 while ( t < r ) {
2753 WORD change;
2754 change = ExtraSymbol(*t,-t[1],nsym,ppsym,&ncoef);
2755 nsym += change;
2756 ppsym += change * 2;
2757 t += 2;
2758 }
2759 }
2760 else {
2761 t += 2;
2762 while ( t < r ) {
2763 *ppdot++ = *t++;
2764 *ppdot++ = *t++;
2765 *ppdot++ = -*t++;
2766 ndot++;
2767 }
2768 }
2769 while ( to < stop ) *to++ = *from++;
2770 r = tt;
2771 }
2772#ifdef WITHFLOAT
2773 else if ( *t == FLOATFUN && TestFloat(t) ) {
2774 k = t[1];
2775 pden[i][1] -= k;
2776 pden[i][FUNHEAD] -= k;
2777 pden[i][FUNHEAD+ARGHEAD] -= k;
2778 UnpackFloat(aux5,t);
2779 if ( withfloat == 0 ) {
2780 mpf_ui_div(aux4,1,aux5);
2781 withfloat = 2;
2782 }
2783 else {
2784 if ( withfloat == 1 ) UnpackFloat(aux4,firstfloat);
2785 mpf_div(aux4,aux4,aux5);
2786 withfloat++;
2787 }
2788 tt = to = t;
2789 while ( r < m ) *to++ = *r++;
2790 *to++ = 1; *to++ = 1; *to++ = 3;
2791 if ( ncoef < 0 ) t[-1] = -t[-1];
2792 m -= k;
2793 stop = m + 3;
2794 r = tt;
2795 }
2796#endif
2797 t = r;
2798 }
2799 if ( pden[i][1] == 4+FUNHEAD+ARGHEAD ) {
2800DropDen:
2801 for ( j = 0; j < nnco; j++ ) {
2802 if ( pden[i] == pnco[j] ) {
2803 --nnco;
2804 while ( j < nnco ) {
2805 pnco[j] = pnco[j+1];
2806 j++;
2807 }
2808 break;
2809 }
2810 }
2811 pden[i--] = pden[--nden];
2812 }
2813 }
2814 }
2815 }
2816/*
2817 #] Easy denominators :
2818 #[ Index Contractions :
2819*/
2820 if ( ndel ) {
2821 t = pdel;
2822 for ( i = 0; i < ndel; i += 2 ) {
2823 if ( t[0] == t[1] ) {
2824 if ( t[0] == EMPTYINDEX ) {}
2825 else if ( *t < AM.OffsetIndex ) {
2826 k = AC.FixIndices[*t];
2827 if ( k < 0 ) { j = -1; k = -k; }
2828 else if ( k > 0 ) j = 1;
2829 else goto NormZero;
2830 goto WithFix;
2831 }
2832 else if ( *t >= AM.DumInd ) {
2833 k = AC.lDefDim;
2834 if ( k ) goto docontract;
2835 }
2836 else if ( *t >= AM.WilInd ) {
2837 k = indices[*t-AM.OffsetIndex-WILDOFFSET].dimension;
2838 if ( k ) goto docontract;
2839 }
2840 else if ( ( k = indices[*t-AM.OffsetIndex].dimension ) != 0 ) {
2841docontract:
2842 if ( k > 0 ) {
2843 j = 1;
2844WithFix: shortnum = k;
2845 ncoef = REDLENG(ncoef);
2846 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),j) )
2847 goto FromNorm;
2848 ncoef = INCLENG(ncoef);
2849 }
2850 else {
2851 WORD change;
2852 change = ExtraSymbol((WORD)(-k),(WORD)1,nsym,ppsym,&ncoef);
2853 nsym += change;
2854 ppsym += change * 2;
2855 }
2856 t[1] = pdel[ndel-1];
2857 t[0] = pdel[ndel-2];
2858HaveCon:
2859 ndel -= 2;
2860 i -= 2;
2861 }
2862 }
2863 else {
2864 if ( *t < AM.OffsetIndex && t[1] < AM.OffsetIndex ) goto NormZero;
2865 j = *t - AM.OffsetIndex;
2866 if ( j >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
2867 || ( *t < AM.WilInd && indices[j].dimension ) ) ) {
2868 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2869 if ( *t == *m ) {
2870 *t = m[1];
2871 *m++ = pdel[ndel-2];
2872 *m = pdel[ndel-1];
2873 goto HaveCon;
2874 }
2875 else if ( *t == m[1] ) {
2876 *t = *m;
2877 *m++ = pdel[ndel-2];
2878 *m = pdel[ndel-1];
2879 goto HaveCon;
2880 }
2881 }
2882 }
2883 j = t[1]-AM.OffsetIndex;
2884 if ( j >= 0 && ( ( t[1] >= AM.DumInd && AC.lDefDim )
2885 || ( t[1] < AM.WilInd && indices[j].dimension ) ) ) {
2886 for ( j = i + 2, m = pdel+j; j < ndel; j += 2, m += 2 ) {
2887 if ( t[1] == *m ) {
2888 t[1] = m[1];
2889 *m++ = pdel[ndel-2];
2890 *m = pdel[ndel-1];
2891 goto HaveCon;
2892 }
2893 else if ( t[1] == m[1] ) {
2894 t[1] = *m;
2895 *m++ = pdel[ndel-2];
2896 *m = pdel[ndel-1];
2897 goto HaveCon;
2898 }
2899 }
2900 }
2901 t += 2;
2902 }
2903 }
2904 if ( ndel > 0 ) {
2905 if ( nvec ) {
2906 t = pdel;
2907 for ( i = 0; i < ndel; i++ ) {
2908 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2909 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2910 r = pvec + 1;
2911 for ( j = 1; j < nvec; j += 2 ) {
2912 if ( *r == *t ) {
2913 if ( i & 1 ) {
2914 *r = t[-1];
2915 *t-- = pdel[--ndel];
2916 i -= 2;
2917 }
2918 else {
2919 *r = t[1];
2920 t[1] = pdel[--ndel];
2921 i--;
2922 }
2923 *t-- = pdel[--ndel];
2924 break;
2925 }
2926 r += 2;
2927 }
2928 }
2929 t++;
2930 }
2931 }
2932 if ( ndel > 0 && ncon ) {
2933 t = pdel;
2934 for ( i = 0; i < ndel; i++ ) {
2935 if ( *t >= AM.OffsetIndex && ( ( *t >= AM.DumInd && AC.lDefDim ) ||
2936 ( *t < AM.WilInd && indices[*t-AM.OffsetIndex].dimension ) ) ) {
2937 for ( j = 0; j < ncon; j++ ) {
2938 if ( *pcon[j] == *t ) {
2939 if ( i & 1 ) {
2940 *pcon[j] = t[-1];
2941 *t-- = pdel[--ndel];
2942 i -= 2;
2943 }
2944 else {
2945 *pcon[j] = t[1];
2946 t[1] = pdel[--ndel];
2947 i--;
2948 }
2949 *t-- = pdel[--ndel];
2950 didcontr++;
2951 r = pcon[j];
2952 for ( j = 0; j < nnco; j++ ) {
2953 m = pnco[j];
2954 if ( r > m && r < m+m[1] ) {
2955 m[2] |= DIRTYSYMFLAG;
2956 break;
2957 }
2958 }
2959 for ( j = 0; j < ncom; j++ ) {
2960 m = pcom[j];
2961 if ( r > m && r < m+m[1] ) {
2962 m[2] |= DIRTYSYMFLAG;
2963 break;
2964 }
2965 }
2966 for ( j = 0; j < neps; j++ ) {
2967 m = peps[j];
2968 if ( r > m && r < m+m[1] ) {
2969 m[2] |= DIRTYSYMFLAG;
2970 break;
2971 }
2972 }
2973 break;
2974 }
2975 }
2976 }
2977 t++;
2978 }
2979 }
2980 }
2981 }
2982 if ( nvec ) {
2983 t = pvec + 1;
2984 for ( i = 3; i < nvec; i += 2 ) {
2985 k = *t - AM.OffsetIndex;
2986 if ( k >= 0 && ( ( *t > AM.DumInd && AC.lDefDim )
2987 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
2988 r = t + 2;
2989 for ( j = i; j < nvec; j += 2 ) {
2990 if ( *r == *t ) { /* Another dotproduct */
2991 *ppdot++ = t[-1];
2992 *ppdot++ = r[-1];
2993 *ppdot++ = 1;
2994 ndot++;
2995 *r-- = pvec[--nvec];
2996 *r = pvec[--nvec];
2997 *t-- = pvec[--nvec];
2998 *t-- = pvec[--nvec];
2999 i -= 2;
3000 break;
3001 }
3002 r += 2;
3003 }
3004 }
3005 t += 2;
3006 }
3007 if ( nvec > 0 && ncon ) {
3008 t = pvec + 1;
3009 for ( i = 1; i < nvec; i += 2 ) {
3010 k = *t - AM.OffsetIndex;
3011 if ( k >= 0 && ( ( *t >= AM.DumInd && AC.lDefDim )
3012 || ( *t < AM.WilInd && indices[k].dimension ) ) ) {
3013 for ( j = 0; j < ncon; j++ ) {
3014 if ( *pcon[j] == *t ) {
3015 *pcon[j] = t[-1];
3016 *t-- = pvec[--nvec];
3017 *t-- = pvec[--nvec];
3018 r = pcon[j];
3019 pcon[j] = pcon[--ncon];
3020 i -= 2;
3021 for ( j = 0; j < nnco; j++ ) {
3022 m = pnco[j];
3023 if ( r > m && r < m+m[1] ) {
3024 m[2] |= DIRTYSYMFLAG;
3025 break;
3026 }
3027 }
3028 for ( j = 0; j < ncom; j++ ) {
3029 m = pcom[j];
3030 if ( r > m && r < m+m[1] ) {
3031 m[2] |= DIRTYSYMFLAG;
3032 break;
3033 }
3034 }
3035 for ( j = 0; j < neps; j++ ) {
3036 m = peps[j];
3037 if ( r > m && r < m+m[1] ) {
3038 m[2] |= DIRTYSYMFLAG;
3039 break;
3040 }
3041 }
3042 break;
3043 }
3044 }
3045 }
3046 t += 2;
3047 }
3048 }
3049 }
3050/*
3051 #] Index Contractions :
3052 #[ NonCommuting Functions :
3053*/
3054 m = fillsetexp;
3055 if ( nnco ) {
3056 for ( i = 0; i < nnco; i++ ) {
3057 t = pnco[i];
3058 if ( ( *t >= (FUNCTION+WILDOFFSET)
3059 && functions[*t-FUNCTION-WILDOFFSET].spec <= 0 )
3060 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3061 && functions[*t-FUNCTION].spec <= 0 ) ) {
3062 DoRevert(t,m);
3063 if ( didcontr ) {
3064 r = t + FUNHEAD;
3065 t += t[1];
3066 while ( r < t ) {
3067 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3068 *r = -SNUMBER;
3069 didcontr--;
3070 pnco[i][2] |= DIRTYSYMFLAG;
3071 }
3072 NEXTARG(r)
3073 }
3074 }
3075 }
3076 }
3077
3078 /* First should come the code for function properties. */
3079
3080 /* First we test for symmetric properties and the DIRTYSYMFLAG */
3081
3082 for ( i = 0; i < nnco; i++ ) {
3083 t = pnco[i];
3084 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) && *t != DOLLAREXPRESSION ) {
3085 l = 0; /* to make the compiler happy */
3086 if ( ( *t >= (FUNCTION+WILDOFFSET)
3087 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3088 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3089 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3090 if ( *t >= (FUNCTION+WILDOFFSET) ) {
3091 *t -= WILDOFFSET;
3092 j = FullSymmetrize(BHEAD t,l);
3093 *t += WILDOFFSET;
3094 }
3095 else j = FullSymmetrize(BHEAD t,l);
3096 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3097 if ( ( j & 2 ) != 0 ) goto NormZero;
3098 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3099 }
3100 }
3101 else t[2] &= ~DIRTYSYMFLAG;
3102 }
3103 }
3104
3105 /* Non commuting functions are then tested for commutation
3106 rules. If needed their order is exchanged. */
3107
3108 k = nnco - 1;
3109 for ( i = 0; i < k; i++ ) {
3110 j = i;
3111 while ( Commute(pnco[j],pnco[j+1]) ) {
3112 t = pnco[j]; pnco[j] = pnco[j+1]; pnco[j+1] = t;
3113 l = j-1;
3114 while ( l >= 0 && Commute(pnco[l],pnco[l+1]) ) {
3115 t = pnco[l]; pnco[l] = pnco[l+1]; pnco[l+1] = t;
3116 l--;
3117 }
3118 if ( ++j >= k ) break;
3119 }
3120 }
3121
3122 /* Finally they are written to output. gamma matrices
3123 are bundled if possible */
3124
3125 for ( i = 0; i < nnco; i++ ) {
3126 t = pnco[i];
3127 if ( *t == IDFUNCTION ) AN.idfunctionflag = 1;
3128 if ( *t >= GAMMA && *t <= GAMMASEVEN ) {
3129 WORD gtype;
3130 to = m;
3131 *m++ = GAMMA;
3132 m++;
3133 FILLFUN(m)
3134 *m++ = stype = t[FUNHEAD]; /* type of string */
3135 j = 0;
3136 nnum = 0;
3137 do {
3138 r = t + t[1];
3139 if ( *t == GAMMAFIVE ) {
3140 gtype = GAMMA5; t += FUNHEAD; goto onegammamatrix; }
3141 else if ( *t == GAMMASIX ) {
3142 gtype = GAMMA6; t += FUNHEAD; goto onegammamatrix; }
3143 else if ( *t == GAMMASEVEN ) {
3144 gtype = GAMMA7; t += FUNHEAD; goto onegammamatrix; }
3145 t += FUNHEAD+1;
3146 while ( t < r ) {
3147 gtype = *t;
3148onegammamatrix:
3149 if ( gtype == GAMMA5 ) {
3150 if ( j == GAMMA1 ) j = GAMMA5;
3151 else if ( j == GAMMA5 ) j = GAMMA1;
3152 else if ( j == GAMMA7 ) ncoef = -ncoef;
3153 if ( nnum & 1 ) ncoef = -ncoef;
3154 }
3155 else if ( gtype == GAMMA6 || gtype == GAMMA7 ) {
3156 if ( nnum & 1 ) {
3157 if ( gtype == GAMMA6 ) gtype = GAMMA7;
3158 else gtype = GAMMA6;
3159 }
3160 if ( j == GAMMA1 ) j = gtype;
3161 else if ( j == GAMMA5 ) {
3162 j = gtype;
3163 if ( j == GAMMA7 ) ncoef = -ncoef;
3164 }
3165 else if ( j != gtype ) goto NormZero;
3166 else {
3167 shortnum = 2;
3168 ncoef = REDLENG(ncoef);
3169 if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(&shortnum),1) ) goto FromNorm;
3170 ncoef = INCLENG(ncoef);
3171 }
3172 }
3173 else {
3174 *m++ = gtype; nnum++;
3175 }
3176 t++;
3177 }
3178
3179 } while ( ( ++i < nnco ) && ( *(t = pnco[i]) >= GAMMA
3180 && *t <= GAMMASEVEN ) && ( t[FUNHEAD] == stype ) );
3181 i--;
3182 if ( j ) {
3183 k = WORDDIF(m,to) - FUNHEAD-1;
3184 r = m;
3185 from = m++;
3186 while ( --k >= 0 ) *from-- = *--r;
3187 *from = j;
3188 }
3189 to[1] = WORDDIF(m,to);
3190 }
3191 else if ( *t < 0 ) {
3192 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3193 FILLFUN3(m)
3194 }
3195 else {
3196 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3197 && *t != REPLACEMENT && *t != DOLLAREXPRESSION
3198 && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3199 k = t[1];
3200 NCOPY(m,t,k);
3201 }
3202 }
3203
3204 }
3205/*
3206 #] NonCommuting Functions :
3207 #[ Commuting Functions :
3208*/
3209 if ( ncom ) {
3210 for ( i = 0; i < ncom; i++ ) {
3211 t = pcom[i];
3212 if ( ( *t >= (FUNCTION+WILDOFFSET)
3213 && functions[*t-FUNCTION-WILDOFFSET].spec <= 0 )
3214 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3215 && functions[*t-FUNCTION].spec <= 0 ) ) {
3216 DoRevert(t,m);
3217 if ( didcontr ) {
3218 r = t + FUNHEAD;
3219 t += t[1];
3220 while ( r < t ) {
3221 if ( *r == -INDEX && r[1] >= 0 && r[1] < AM.OffsetIndex ) {
3222 *r = -SNUMBER;
3223 didcontr--;
3224 pcom[i][2] |= DIRTYSYMFLAG;
3225 }
3226 NEXTARG(r)
3227 }
3228 }
3229 }
3230 }
3231
3232 /* Now we test for symmetric properties and the DIRTYSYMFLAG */
3233
3234 for ( i = 0; i < ncom; i++ ) {
3235 t = pcom[i];
3236 if ( *t > 0 && ( t[2] & DIRTYSYMFLAG ) ) {
3237 l = 0; /* to make the compiler happy */
3238 if ( ( *t >= (FUNCTION+WILDOFFSET)
3239 && ( l = functions[*t-FUNCTION-WILDOFFSET].symmetric ) > 0 )
3240 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3241 && ( l = functions[*t-FUNCTION].symmetric ) > 0 ) ) {
3242 if ( *t >= (FUNCTION+WILDOFFSET) ) {
3243 *t -= WILDOFFSET;
3244 j = FullSymmetrize(BHEAD t,l);
3245 *t += WILDOFFSET;
3246 }
3247 else j = FullSymmetrize(BHEAD t,l);
3248 if ( (l & ~REVERSEORDER) == ANTISYMMETRIC ) {
3249 if ( ( j & 2 ) != 0 ) goto NormZero;
3250 if ( ( j & 1 ) != 0 ) ncoef = -ncoef;
3251 }
3252 }
3253 else t[2] &= ~DIRTYSYMFLAG;
3254 }
3255 }
3256/*
3257 Sort the functions
3258 From a purists point of view this can be improved.
3259 There are slow and fast arguments and no conversions are
3260 taken into account here.
3261*/
3262 for ( i = 1; i < ncom; i++ ) {
3263 for ( j = i; j > 0; j-- ) {
3264 WORD jj,kk;
3265 jj = j-1;
3266 t = pcom[jj];
3267 r = pcom[j];
3268 if ( *t < 0 ) {
3269 if ( *r < 0 ) { if ( *t >= *r ) goto NextI; }
3270 else { if ( -*t <= *r ) goto NextI; }
3271 goto jexch;
3272 }
3273 else if ( *r < 0 ) {
3274 if ( *t < -*r ) goto NextI;
3275 goto jexch;
3276 }
3277 else if ( *t != *r ) {
3278 if ( *t < *r ) goto NextI;
3279jexch: t = pcom[j]; pcom[j] = pcom[jj]; pcom[jj] = t;
3280 continue;
3281 }
3282 if ( AC.properorderflag ) {
3283 if ( ( *t >= (FUNCTION+WILDOFFSET)
3284 && functions[*t-FUNCTION-WILDOFFSET].spec >= TENSORFUNCTION )
3285 || ( *t >= FUNCTION && *t < (FUNCTION + WILDOFFSET)
3286 && functions[*t-FUNCTION].spec >= TENSORFUNCTION ) ) {}
3287 else {
3288 WORD *s1, *s2, *ss1, *ss2;
3289 s1 = t+FUNHEAD; s2 = r+FUNHEAD;
3290 ss1 = t + t[1]; ss2 = r + r[1];
3291 while ( s1 < ss1 && s2 < ss2 ) {
3292 k = CompArg(s1,s2);
3293 if ( k > 0 ) goto jexch;
3294 if ( k < 0 ) goto NextI;
3295 NEXTARG(s1)
3296 NEXTARG(s2)
3297 }
3298 if ( s1 < ss1 ) goto jexch;
3299 goto NextI;
3300 }
3301 k = t[1] - FUNHEAD;
3302 kk = r[1] - FUNHEAD;
3303 t += FUNHEAD;
3304 r += FUNHEAD;
3305 while ( k > 0 && kk > 0 ) {
3306 if ( *t < *r ) goto NextI;
3307 else if ( *t++ > *r++ ) goto jexch;
3308 k--; kk--;
3309 }
3310 if ( k > 0 ) goto jexch;
3311 goto NextI;
3312 }
3313 else
3314 {
3315 k = t[1] - FUNHEAD;
3316 kk = r[1] - FUNHEAD;
3317 t += FUNHEAD;
3318 r += FUNHEAD;
3319 while ( k > 0 && kk > 0 ) {
3320 if ( *t < *r ) goto NextI;
3321 else if ( *t++ > *r++ ) goto jexch;
3322 k--; kk--;
3323 }
3324 if ( k > 0 ) goto jexch;
3325 goto NextI;
3326 }
3327 }
3328NextI:;
3329 }
3330 for ( i = 0; i < ncom; i++ ) {
3331 t = pcom[i];
3332 if ( *t == THETA || *t == THETA2 ) {
3333 if ( ( k = DoTheta(BHEAD t) ) == 0 ) goto NormZero;
3334 else if ( k < 0 ) {
3335 k = t[1];
3336 NCOPY(m,t,k);
3337 }
3338 }
3339 else if ( *t == DELTA2 || *t == DELTAP ) {
3340 if ( ( k = DoDelta(t) ) == 0 ) goto NormZero;
3341 else if ( k < 0 ) {
3342 k = t[1];
3343 NCOPY(m,t,k);
3344 }
3345 }
3346 else if ( *t == AR.PolyFunInv && AR.PolyFunType == 2 ) {
3347/*
3348 If there are two arguments, exchange them, change the
3349 name of the function and go to dealing with PolyRatFun.
3350*/
3351 WORD *mm, *tt = t, numt = 0;
3352 tt += FUNHEAD;
3353 while ( tt < t+t[1] ) { numt++; NEXTARG(tt) }
3354 if ( numt == 2 ) {
3355 tt = t; mm = m; k = t[1];
3356 NCOPY(mm,tt,k)
3357 mm = m+FUNHEAD;
3358 NEXTARG(mm);
3359 tt = t+FUNHEAD;
3360 if ( *mm < 0 ) {
3361 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3362 else { *tt++ = *mm++; *tt++ = *mm++; }
3363 }
3364 else {
3365 k = *mm; NCOPY(tt,mm,k)
3366 }
3367 mm = m+FUNHEAD;
3368 if ( *mm < 0 ) {
3369 if ( *mm <= -FUNCTION ) { *tt++ = *mm++; }
3370 else { *tt++ = *mm++; *tt++ = *mm++; }
3371 }
3372 else {
3373 k = *mm; NCOPY(tt,mm,k)
3374 }
3375 *t = AR.PolyFun;
3376 t[2] |= MUSTCLEANPRF;
3377 goto regularratfun;
3378 }
3379 }
3380 else if ( *t == AR.PolyFun ) {
3381 if ( AR.PolyFunType == 1 ) { /* Regular PolyFun with one argument */
3382 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3383 t[1] == FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) goto NormZero;
3384 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) {
3385 if ( AN.PolyNormFlag == 0 ) {
3386 AN.PolyNormFlag = 1;
3387 AN.PolyFunTodo = 0;
3388 }
3389 }
3390 k = t[1];
3391 NCOPY(m,t,k);
3392 }
3393 else if ( AR.PolyFunType == 2 ) {
3394/*
3395 PolyRatFun.
3396 Regular type: Two arguments
3397 Power expanded: One argument. Here to be treated as
3398 AR.PolyFunType == 1, but with power cutoff.
3399*/
3400regularratfun:;
3401/*
3402 First check for zeroes.
3403*/
3404 if ( t[FUNHEAD+1] == 0 && AR.Eside != LHSIDE &&
3405 t[1] > FUNHEAD + 2 && t[FUNHEAD] == -SNUMBER ) {
3406 u = t + FUNHEAD + 2;
3407 if ( *u < 0 ) {
3408 if ( *u <= -FUNCTION ) {}
3409 else if ( t[1] == FUNHEAD+4 && t[FUNHEAD+2] == -SNUMBER
3410 && t[FUNHEAD+3] == 0 ) goto NormPRF;
3411 else if ( t[1] == FUNHEAD+4 ) goto NormZero;
3412 }
3413 else if ( t[1] == *u+FUNHEAD+2 ) goto NormZero;
3414 }
3415 else {
3416 u = t+FUNHEAD; NEXTARG(u);
3417 if ( *u == -SNUMBER && u[1] == 0 ) goto NormInf;
3418 }
3419 if ( i > 0 && pcom[i-1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3420 else if ( i < ncom-1 && pcom[i+1][0] == AR.PolyFun ) AN.PolyNormFlag = 1;
3421 k = t[1];
3422 if ( AN.PolyNormFlag ) {
3423 if ( AR.PolyFunExp == 0 ) {
3424 AN.PolyFunTodo = 0;
3425 NCOPY(m,t,k);
3426 }
3427 else if ( AR.PolyFunExp == 1 ) { /* get highest divergence */
3428 if ( PolyFunMode == 0 ) {
3429 NCOPY(m,t,k);
3430 AN.PolyFunTodo = 1;
3431 }
3432 else {
3433 WORD *mmm = m;
3434 NCOPY(m,t,k);
3435 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3436 goto FromNorm;
3437 m = mmm+mmm[1];
3438 }
3439 }
3440 else {
3441 if ( PolyFunMode == 0 ) {
3442 NCOPY(m,t,k);
3443 AN.PolyFunTodo = 1;
3444 }
3445 else {
3446 WORD *mmm = m;
3447 NCOPY(m,t,k);
3448 if ( ExpandRat(BHEAD mmm) != 0 )
3449 goto FromNorm;
3450 m = mmm+mmm[1];
3451 }
3452 }
3453 }
3454 else {
3455 if ( AR.PolyFunExp == 0 ) {
3456 AN.PolyFunTodo = 0;
3457 NCOPY(m,t,k);
3458 }
3459 else if ( AR.PolyFunExp == 1 ) { /* get highest divergence */
3460 WORD *mmm = m;
3461 NCOPY(m,t,k);
3462 if ( TreatPolyRatFun(BHEAD mmm) != 0 )
3463 goto FromNorm;
3464 m = mmm+mmm[1];
3465 }
3466 else {
3467 WORD *mmm = m;
3468 NCOPY(m,t,k);
3469 if ( ExpandRat(BHEAD mmm) != 0 )
3470 goto FromNorm;
3471 m = mmm+mmm[1];
3472 }
3473 }
3474 }
3475 }
3476 else if ( *t > 0 ) {
3477 if ( ( t[2] & DIRTYFLAG ) == DIRTYFLAG
3478 && *t != REPLACEMENT && TestFunFlag(BHEAD t) ) ReplaceVeto = 1;
3479 k = t[1];
3480 NCOPY(m,t,k);
3481 }
3482 else {
3483 *m++ = -*t; *m++ = FUNHEAD; *m++ = 0;
3484 FILLFUN3(m)
3485 }
3486 }
3487 }
3488/*
3489 #] Commuting Functions :
3490 #[ Track Replace_ :
3491*/
3492 if ( ReplaceVeto < 0 ) {
3493/*
3494 We found one (or more) replace_ functions and all other
3495 functions are 'clean' (no dirty flag).
3496 Now we check whether one of these functions can be used.
3497 Thus far the functions go from fillsetexp to m.
3498 Somewhere in there there are -ReplaceVeto occurrences of REPLACEMENT.
3499 Hunt for the first one that fits the bill.
3500 Note that replace_ is a commuting function.
3501*/
3502 WORD *ma = fillsetexp, *mb, *mc;
3503 while ( ma < m ) {
3504 mb = ma + ma[1];
3505 if ( *ma != REPLACEMENT ) {
3506 ma = mb;
3507 continue;
3508 }
3509 if ( *ma == REPLACEMENT && ReplaceType == -1 ) {
3510 mc = ma;
3511 ReplaceType = 0;
3512 if ( AN.RSsize < 2*ma[1]+SUBEXPSIZE ) {
3513 if ( AN.ReplaceScrat ) M_free(AN.ReplaceScrat,"AN.ReplaceScrat");
3514 AN.RSsize = 2*ma[1]+SUBEXPSIZE+40;
3515 AN.ReplaceScrat = (WORD *)Malloc1((AN.RSsize+1)*sizeof(WORD),"AN.ReplaceScrat");
3516 }
3517 ma += FUNHEAD;
3518 ReplaceSub = AN.ReplaceScrat;
3519 ReplaceSub += SUBEXPSIZE;
3520 while ( ma < mb ) {
3521 if ( *ma > 0 ) goto NoRep;
3522 if ( *ma <= -FUNCTION ) {
3523 *ReplaceSub++ = FUNTOFUN;
3524 *ReplaceSub++ = 4;
3525 *ReplaceSub++ = -*ma++;
3526 if ( *ma > -FUNCTION ) goto NoRep;
3527 *ReplaceSub++ = -*ma++;
3528 }
3529 else if ( ma+4 > mb ) goto NoRep;
3530 else {
3531 if ( *ma == -SYMBOL ) {
3532 if ( ma[2] == -SYMBOL && ma+4 <= mb )
3533 *ReplaceSub++ = SYMTOSYM;
3534 else if ( ma[2] == -SNUMBER && ma+4 <= mb ) {
3535 *ReplaceSub++ = SYMTONUM;
3536 if ( ReplaceType == 0 ) {
3537 oldtoprhs = C->numrhs;
3538 oldcpointer = C->Pointer - C->Buffer;
3539 }
3540 ReplaceType = 1;
3541 }
3542 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3543 *ReplaceSub++ = SYMTONUM;
3544 *ReplaceSub++ = 4;
3545 *ReplaceSub++ = ma[1];
3546 *ReplaceSub++ = 0;
3547 ma += 2+ARGHEAD;
3548 continue;
3549 }
3550/*
3551 Next is the subexpression. We have to test that
3552 it isn't vector-like or index-like
3553*/
3554 else if ( ma[2] > 0 ) {
3555 WORD *sstop, *ttstop, n;
3556 ss = ma+2;
3557 sstop = ss + *ss;
3558 ss += ARGHEAD;
3559 while ( ss < sstop ) {
3560 tt = ss + *ss;
3561 ttstop = tt - ABS(tt[-1]);
3562 ss++;
3563 while ( ss < ttstop ) {
3564 if ( *ss == INDEX ) goto NoRep;
3565 ss += ss[1];
3566 }
3567 ss = tt;
3568 }
3569 subtype = SYMTOSUB;
3570 if ( ReplaceType == 0 ) {
3571 oldtoprhs = C->numrhs;
3572 oldcpointer = C->Pointer - C->Buffer;
3573 }
3574 ReplaceType = 1;
3575 ss = AddRHS(AT.ebufnum,1);
3576 tt = ma+2;
3577 n = *tt - ARGHEAD;
3578 tt += ARGHEAD;
3579 while ( (ss + n + 10) > C->Top ) ss = DoubleCbuffer(AT.ebufnum,ss,14);
3580 while ( --n >= 0 ) *ss++ = *tt++;
3581 *ss++ = 0;
3582 C->rhs[C->numrhs+1] = ss;
3583 C->Pointer = ss;
3584 *ReplaceSub++ = subtype;
3585 *ReplaceSub++ = 4;
3586 *ReplaceSub++ = ma[1];
3587 *ReplaceSub++ = C->numrhs;
3588 ma += 2 + ma[2];
3589 continue;
3590 }
3591 else goto NoRep;
3592 }
3593 else if ( ( *ma == -VECTOR || *ma == -MINVECTOR ) && ma+4 <= mb ) {
3594 if ( ma[2] == -VECTOR ) {
3595 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOVEC;
3596 else *ReplaceSub++ = VECTOMIN;
3597 }
3598 else if ( ma[2] == -MINVECTOR ) {
3599 if ( *ma == -VECTOR ) *ReplaceSub++ = VECTOMIN;
3600 else *ReplaceSub++ = VECTOVEC;
3601 }
3602/*
3603 Next is a vector-like subexpression
3604 Search for vector nature first
3605*/
3606 else if ( ma[2] > 0 ) {
3607 WORD *sstop, *ttstop, *w, *mm, n, count;
3608 WORD *v1, *v2 = 0;
3609 if ( *ma == -MINVECTOR ) {
3610 ss = ma+2;
3611 sstop = ss + *ss;
3612 ss += ARGHEAD;
3613 while ( ss < sstop ) {
3614 ss += *ss;
3615 ss[-1] = -ss[-1];
3616 }
3617 *ma = -VECTOR;
3618 }
3619 ss = ma+2;
3620 sstop = ss + *ss;
3621 ss += ARGHEAD;
3622 while ( ss < sstop ) {
3623 tt = ss + *ss;
3624 ttstop = tt - ABS(tt[-1]);
3625 ss++;
3626 count = 0;
3627 while ( ss < ttstop ) {
3628 if ( *ss == INDEX ) {
3629 n = ss[1] - 2; ss += 2;
3630 while ( --n >= 0 ) {
3631 if ( *ss < MINSPEC ) count++;
3632 ss++;
3633 }
3634 }
3635 else ss += ss[1];
3636 }
3637 if ( count != 1 ) goto NoRep;
3638 ss = tt;
3639 }
3640 subtype = VECTOSUB;
3641 if ( ReplaceType == 0 ) {
3642 oldtoprhs = C->numrhs;
3643 oldcpointer = C->Pointer - C->Buffer;
3644 }
3645 ReplaceType = 1;
3646 mm = AddRHS(AT.ebufnum,1);
3647 *ReplaceSub++ = subtype;
3648 *ReplaceSub++ = 4;
3649 *ReplaceSub++ = ma[1];
3650 *ReplaceSub++ = C->numrhs;
3651 w = ma+2;
3652 n = *w - ARGHEAD;
3653 w += ARGHEAD;
3654 while ( (mm + n + 10) > C->Top )
3655 mm = DoubleCbuffer(AT.ebufnum,mm,15);
3656 while ( --n >= 0 ) *mm++ = *w++;
3657 *mm++ = 0;
3658 C->rhs[C->numrhs+1] = mm;
3659 C->Pointer = mm;
3660 mm = AddRHS(AT.ebufnum,1);
3661 w = ma+2;
3662 n = *w - ARGHEAD;
3663 w += ARGHEAD;
3664 while ( (mm + n + 13) > C->Top )
3665 mm = DoubleCbuffer(AT.ebufnum,mm,16);
3666 sstop = w + n;
3667 while ( w < sstop ) {
3668 tt = w + *w; ttstop = tt - ABS(tt[-1]);
3669 ss = mm; mm++; w++;
3670 while ( w < ttstop ) { /* Subterms */
3671 if ( *w != INDEX ) {
3672 n = w[1];
3673 NCOPY(mm,w,n);
3674 }
3675 else {
3676 v1 = mm;
3677 *mm++ = *w++;
3678 *mm++ = n = *w++;
3679 n -= 2;
3680 while ( --n >= 0 ) {
3681 if ( *w >= MINSPEC ) *mm++ = *w++;
3682 else v2 = w++;
3683 }
3684 n = WORDDIF(mm,v1);
3685 if ( n != v1[1] ) {
3686 if ( n <= 2 ) mm -= 2;
3687 else v1[1] = n;
3688 *mm++ = VECTOR;
3689 *mm++ = 4;
3690 *mm++ = *v2;
3691 *mm++ = FUNNYVEC;
3692 }
3693 }
3694 }
3695 while ( w < tt ) *mm++ = *w++;
3696 *ss = WORDDIF(mm,ss);
3697 }
3698 *mm++ = 0;
3699 C->rhs[C->numrhs+1] = mm;
3700 C->Pointer = mm;
3701 if ( mm > C->Top ) {
3702 MLOCK(ErrorMessageLock);
3703 MesPrint("Internal error in Normalize with extra compiler buffer");
3704 MUNLOCK(ErrorMessageLock);
3705 Terminate(-1);
3706 }
3707 ma += 2 + ma[2];
3708 continue;
3709 }
3710 else goto NoRep;
3711 }
3712 else if ( *ma == -INDEX ) {
3713 if ( ( ma[2] == -INDEX || ma[2] == -VECTOR )
3714 && ma+4 <= mb )
3715 *ReplaceSub++ = INDTOIND;
3716 else if ( ma[1] >= AM.OffsetIndex ) {
3717 if ( ma[2] == -SNUMBER && ma+4 <= mb
3718 && ma[3] >= 0 && ma[3] < AM.OffsetIndex )
3719 *ReplaceSub++ = INDTOIND;
3720 else if ( ma[2] == ARGHEAD && ma+2+ARGHEAD <= mb ) {
3721 *ReplaceSub++ = INDTOIND;
3722 *ReplaceSub++ = 4;
3723 *ReplaceSub++ = ma[1];
3724 *ReplaceSub++ = 0;
3725 ma += 2+ARGHEAD;
3726 continue;
3727 }
3728 else goto NoRep;
3729 }
3730 else goto NoRep;
3731 }
3732 else goto NoRep;
3733 *ReplaceSub++ = 4;
3734 *ReplaceSub++ = ma[1];
3735 *ReplaceSub++ = ma[3];
3736 ma += 4;
3737 }
3738
3739 }
3740 AN.ReplaceScrat[1] = ReplaceSub-AN.ReplaceScrat;
3741/*
3742 Success. This means that we have to remove the replace_
3743 from the functions. It starts at mc and end at mb.
3744*/
3745 while ( mb < m ) *mc++ = *mb++;
3746 m = mc;
3747 break;
3748NoRep:
3749 if ( ReplaceType > 0 ) {
3750 C->numrhs = oldtoprhs;
3751 C->Pointer = C->Buffer + oldcpointer;
3752 }
3753 ReplaceType = -1;
3754 if ( ++ReplaceVeto >= 0 ) break;
3755 }
3756 ma = mb;
3757 }
3758 }
3759/*
3760 #] Track Replace_ :
3761 #[ LeviCivita tensors :
3762*/
3763 if ( neps ) {
3764 to = m;
3765 for ( i = 0; i < neps; i++ ) { /* Put the indices in order */
3766 t = peps[i];
3767 if ( ( t[2] & DIRTYSYMFLAG ) != DIRTYSYMFLAG ) continue;
3768 t[2] &= ~DIRTYSYMFLAG;
3769 if ( AR.Eside == LHSIDE || AR.Eside == LHSIDEX ) {
3770 /* Potential problems with FUNNYWILD */
3771/*
3772 First make sure all FUNNIES are at the end.
3773 Then sort separately
3774*/
3775 r = t + FUNHEAD;
3776 m = tt = t + t[1];
3777 while ( r < m ) {
3778 if ( *r != FUNNYWILD ) { r++; continue; }
3779 k = r[1]; u = r + 2;
3780 while ( u < tt ) {
3781 u[-2] = *u;
3782 if ( *u != FUNNYWILD ) ncoef = -ncoef;
3783 u++;
3784 }
3785 tt[-2] = FUNNYWILD; tt[-1] = k; m -= 2;
3786 }
3787 t += FUNHEAD;
3788 do {
3789 for ( r = t + 1; r < m; r++ ) {
3790 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3791 else if ( *r == *t ) goto NormZero;
3792 }
3793 t++;
3794 } while ( t < m );
3795 do {
3796 for ( r = t + 2; r < tt; r += 2 ) {
3797 if ( r[1] < t[1] ) {
3798 k = r[1]; r[1] = t[1]; t[1] = k; ncoef = -ncoef; }
3799 else if ( r[1] == t[1] ) goto NormZero;
3800 }
3801 t += 2;
3802 } while ( t < tt );
3803 }
3804 else {
3805 m = t + t[1];
3806 t += FUNHEAD;
3807 do {
3808 for ( r = t + 1; r < m; r++ ) {
3809 if ( *r < *t ) { k = *r; *r = *t; *t = k; ncoef = -ncoef; }
3810 else if ( *r == *t ) goto NormZero;
3811 }
3812 t++;
3813 } while ( t < m );
3814 }
3815 }
3816
3817 /* Sort the tensors */
3818
3819 for ( i = 0; i < (neps-1); i++ ) {
3820 t = peps[i];
3821 for ( j = i+1; j < neps; j++ ) {
3822 r = peps[j];
3823 if ( t[1] > r[1] ) {
3824 peps[i] = m = r; peps[j] = r = t; t = m;
3825 }
3826 else if ( t[1] == r[1] ) {
3827 k = t[1] - FUNHEAD;
3828 m = t + FUNHEAD;
3829 r += FUNHEAD;
3830 do {
3831 if ( *r < *m ) {
3832 m = peps[j]; peps[j] = t; peps[i] = t = m;
3833 break;
3834 }
3835 else if ( *r++ > *m++ ) break;
3836 } while ( --k > 0 );
3837 }
3838 }
3839 }
3840 m = to;
3841 for ( i = 0; i < neps; i++ ) {
3842 t = peps[i];
3843 k = t[1];
3844 NCOPY(m,t,k);
3845 }
3846 }
3847/*
3848 #] LeviCivita tensors :
3849 #[ Delta :
3850*/
3851 if ( ndel ) {
3852 r = t = pdel;
3853 for ( i = 0; i < ndel; i += 2, r += 2 ) {
3854 if ( r[1] < r[0] ) { k = *r; *r = r[1]; r[1] = k; }
3855 }
3856 for ( i = 2; i < ndel; i += 2, t += 2 ) {
3857 r = t + 2;
3858 for ( j = i; j < ndel; j += 2 ) {
3859 if ( *r > *t ) { r += 2; }
3860 else if ( *r < *t ) {
3861 k = *r; *r++ = *t; *t++ = k;
3862 k = *r; *r++ = *t; *t-- = k;
3863 }
3864 else {
3865 if ( *++r < t[1] ) {
3866 k = *r; *r = t[1]; t[1] = k;
3867 }
3868 r++;
3869 }
3870 }
3871 }
3872 t = pdel;
3873 *m++ = DELTA;
3874 *m++ = ndel + 2;
3875 i = ndel;
3876 NCOPY(m,t,i);
3877 }
3878/*
3879 #] Delta :
3880 #[ Loose Vectors/Indices :
3881*/
3882 if ( nind ) {
3883 t = pind;
3884 for ( i = 0; i < nind; i++ ) {
3885 r = t + 1;
3886 for ( j = i+1; j < nind; j++ ) {
3887 if ( *r < *t ) {
3888 k = *r; *r = *t; *t = k;
3889 }
3890 r++;
3891 }
3892 t++;
3893 }
3894 t = pind;
3895 *m++ = INDEX;
3896 *m++ = nind + 2;
3897 i = nind;
3898 NCOPY(m,t,i);
3899 }
3900/*
3901 #] Loose Vectors/Indices :
3902 #[ Vectors :
3903*/
3904 if ( nvec ) {
3905 t = pvec;
3906 for ( i = 2; i < nvec; i += 2 ) {
3907 r = t + 2;
3908 for ( j = i; j < nvec; j += 2 ) {
3909 if ( *r == *t ) {
3910 if ( *++r < t[1] ) {
3911 k = *r; *r = t[1]; t[1] = k;
3912 }
3913 r++;
3914 }
3915 else if ( *r < *t ) {
3916 k = *r; *r++ = *t; *t++ = k;
3917 k = *r; *r++ = *t; *t-- = k;
3918 }
3919 else { r += 2; }
3920 }
3921 t += 2;
3922 }
3923 t = pvec;
3924 *m++ = VECTOR;
3925 *m++ = nvec + 2;
3926 i = nvec;
3927 NCOPY(m,t,i);
3928 }
3929/*
3930 #] Vectors :
3931 #[ Dotproducts :
3932*/
3933 if ( ndot ) {
3934 to = m;
3935 m = t = pdot;
3936 i = ndot;
3937 while ( --i >= 0 ) {
3938 if ( *t > t[1] ) { j = *t; *t = t[1]; t[1] = j; }
3939 t += 3;
3940 }
3941 t = m;
3942 ndot *= 3;
3943 m += ndot;
3944 while ( t < (m-3) ) {
3945 r = t + 3;
3946 do {
3947 if ( *r == *t ) {
3948 if ( *++r == *++t ) {
3949 r++;
3950 if ( ( *r < MAXPOWER && t[1] < MAXPOWER )
3951 || ( *r > -MAXPOWER && t[1] > -MAXPOWER ) ) {
3952 t++;
3953 *t += *r;
3954 if ( *t > MAXPOWER || *t < -MAXPOWER ) {
3955 MLOCK(ErrorMessageLock);
3956 MesPrint("Exponent of dotproduct out of range: %d",*t);
3957 MUNLOCK(ErrorMessageLock);
3958 goto NormMin;
3959 }
3960 ndot -= 3;
3961 *r-- = *--m;
3962 *r-- = *--m;
3963 *r = *--m;
3964 if ( !*t ) {
3965 ndot -= 3;
3966 *t-- = *--m;
3967 *t-- = *--m;
3968 *t = *--m;
3969 t -= 3;
3970 break;
3971 }
3972 }
3973 else if ( *r < *++t ) {
3974 k = *r; *r++ = *t; *t = k;
3975 }
3976 else r++;
3977 t -= 2;
3978 }
3979 else if ( *r < *t ) {
3980 k = *r; *r++ = *t; *t++ = k;
3981 k = *r; *r++ = *t; *t = k;
3982 t -= 2;
3983 }
3984 else { r += 2; t--; }
3985 }
3986 else if ( *r < *t ) {
3987 k = *r; *r++ = *t; *t++ = k;
3988 k = *r; *r++ = *t; *t++ = k;
3989 k = *r; *r++ = *t; *t = k;
3990 t -= 2;
3991 }
3992 else { r += 3; }
3993 } while ( r < m );
3994 t += 3;
3995 }
3996 m = to;
3997 t = pdot;
3998 if ( ( i = ndot ) > 0 ) {
3999 *m++ = DOTPRODUCT;
4000 *m++ = i + 2;
4001 NCOPY(m,t,i);
4002 }
4003 }
4004/*
4005 #] Dotproducts :
4006 #[ Symbols :
4007*/
4008 if ( nsym ) {
4009 nsym <<= 1;
4010 t = psym;
4011 *m++ = SYMBOL;
4012 r = m;
4013 *m++ = ( i = nsym ) + 2;
4014 if ( i ) { do {
4015 if ( !*t ) {
4016 if ( t[1] < (2*MAXPOWER) ) { /* powers of i */
4017 if ( t[1] & 1 ) { *m++ = 0; *m++ = 1; }
4018 else *r -= 2;
4019 if ( *++t & 2 ) ncoef = -ncoef;
4020 t++;
4021 }
4022 }
4023 else if ( *t <= NumSymbols && *t > -2*MAXPOWER ) { /* Put powers in range */
4024 if ( ( ( ( t[1] > symbols[*t].maxpower ) && ( symbols[*t].maxpower < MAXPOWER ) ) ||
4025 ( ( t[1] < symbols[*t].minpower ) && ( symbols[*t].minpower > -MAXPOWER ) ) ) &&
4026 ( t[1] < 2*MAXPOWER ) && ( t[1] > -2*MAXPOWER ) ) {
4027 if ( i <= 2 || t[2] != *t ) goto NormZero;
4028 }
4029 if ( AN.ncmod == 1 && ( AC.modmode & ALSOPOWERS ) != 0 ) {
4030 if ( AC.cmod[0] == 1 ) t[1] = 0;
4031 else if ( t[1] >= 0 ) t[1] = 1 + (t[1]-1)%(AC.cmod[0]-1);
4032 else {
4033 t[1] = -1 - (-t[1]-1)%(AC.cmod[0]-1);
4034 if ( t[1] < 0 ) t[1] += (AC.cmod[0]-1);
4035 }
4036 }
4037 if ( ( t[1] < (2*MAXPOWER) && t[1] >= MAXPOWER )
4038 || ( t[1] > -(2*MAXPOWER) && t[1] <= -MAXPOWER ) ) {
4039 MLOCK(ErrorMessageLock);
4040 MesPrint("Exponent out of range: %d",t[1]);
4041 MUNLOCK(ErrorMessageLock);
4042 goto NormMin;
4043 }
4044 if ( AT.TrimPower && AR.PolyFunVar == *t && t[1] > AR.PolyFunPow ) {
4045 goto NormZero;
4046 }
4047 else if ( t[1] ) {
4048 *m++ = *t++;
4049 *m++ = *t++;
4050 }
4051 else { *r -= 2; t += 2; }
4052 }
4053 else {
4054 *m++ = *t++; *m++ = *t++;
4055 }
4056 } while ( (i-=2) > 0 ); }
4057 if ( *r <= 2 ) m = r-1;
4058 }
4059/*
4060 #] Symbols :
4061 #[ float_ :
4062
4063 Here we treat float_ functions and combined them with the regular
4064 coefficient.
4065*/
4066
4067#ifdef WITHFLOAT
4068 if ( withfloat ) {
4069 WORD floatsign = 3;
4070/*
4071 First check whether the coefficient is already 1/1
4072*/
4073 if ( ABS(ncoef) == 3 && n_coef[0] == 1 && n_coef[1] == 1 ) {
4074 if ( withfloat == 1 ) {
4075 t = firstfloat;
4076/*
4077 We only interfere by making the sign positive.
4078 The float_ function has already been tested.
4079*/
4080 if ( t[FUNHEAD+3] < 0 ) {
4081 t[FUNHEAD+3] = -t[FUNHEAD+3];
4082 floatsign = -floatsign;
4083 }
4084 if ( ncoef < 0 ) floatsign = -floatsign;
4085/*
4086 Copy to output;
4087*/
4088 AT.FloatPos = m-termout;
4089 i = t[1]; NCOPY(m,t,i)
4090 }
4091 else {
4092 PackFloat(m,aux4);
4093 if ( m[FUNHEAD+3] < 0 ) {
4094 m[FUNHEAD+3] = -m[FUNHEAD+3];
4095 floatsign = -floatsign;
4096 }
4097 if ( ncoef < 0 ) floatsign = -floatsign;
4098 AT.FloatPos = m-termout;
4099 m += m[1];
4100 }
4101 }
4102 else {
4103 if ( withfloat == 1 ) UnpackFloat(aux4,firstfloat);
4104 RatToFloat(aux5,(UWORD *)n_coef,ncoef);
4105 mpf_mul(aux4,aux4,aux5);
4106 PackFloat(m,aux4);
4107 AT.FloatPos = m-termout;
4108 if ( m[FUNHEAD+3] < 0 ) {
4109 m[FUNHEAD+3] = -m[FUNHEAD+3];
4110 floatsign = -floatsign;
4111 }
4112 m += m[1];
4113 }
4114 n_coef[0] = 1; n_coef[1] = 1; ncoef = floatsign;
4115 }
4116 else AT.FloatPos = 0;
4117#endif
4118
4119/*
4120 #] float_ :
4121 #[ Do Replace_ :
4122*/
4123 stop = (WORD *)(((UBYTE *)(termout)) + AM.MaxTer);
4124 i = ABS(ncoef);
4125 if ( ( m + i ) > stop ) {
4126 MLOCK(ErrorMessageLock);
4127 MesPrint("Term too complex during normalization");
4128 MUNLOCK(ErrorMessageLock);
4129 goto NormMin;
4130 }
4131 if ( AT.SS == AT.S0 ) {
4132 if ( ( m + i - termout ) > AT.SS->verbMaxTermSize ) {
4133 AT.SS->verbMaxTermSize = m+i-termout;
4134 }
4135 }
4136 if ( ReplaceType >= 0 ) {
4137 t = n_coef;
4138 i--;
4139 NCOPY(m,t,i);
4140 *m++ = ncoef;
4141 t = termout;
4142 *t = WORDDIF(m,t);
4143 if ( ReplaceType == 0 ) {
4144 AT.WorkPointer = termout+*termout;
4145 WildFill(BHEAD term,termout,AN.ReplaceScrat);
4146 termout = term + *term;
4147 }
4148 else {
4149 AT.WorkPointer = r = termout + *termout;
4150 WildFill(BHEAD r,termout,AN.ReplaceScrat);
4151 i = *r; m = term;
4152 NCOPY(m,r,i);
4153 termout = m;
4154
4155
4156 r = m = term;
4157 r += *term; r -= ABS(r[-1]);
4158 m++;
4159 while ( m < r ) {
4160 if ( *m >= FUNCTION && m[1] > FUNHEAD &&
4161 functions[*m-FUNCTION].spec != TENSORFUNCTION )
4162 m[2] |= DIRTYFLAG;
4163 m += m[1];
4164 }
4165 }
4166/*
4167 The next 'reset' cannot be done. We still need the expression
4168 in the buffer. Note though that this may cause a runaway pointer
4169 if we are not very careful.
4170
4171 C->numrhs = oldtoprhs;
4172 C->Pointer = C->Buffer + oldcpointer;
4173*/
4174 AT.NormDepth--;
4175 TermFree(n_llnum,"n_llnum");
4176 TermFree(n_coef,"NormCoef");
4177 return(1);
4178 }
4179 else {
4180 t = termout;
4181 k = WORDDIF(m,t);
4182 *t = k + i;
4183 m = term;
4184 NCOPY(m,t,k);
4185 i--;
4186 t = n_coef;
4187 NCOPY(m,t,i);
4188 *m++ = ncoef;
4189 }
4190/*
4191 #] Do Replace_ :
4192 #[ Errors and Finish :
4193*/
4194RegEnd:
4195 if ( termout < term + *term && termout >= term ) AT.WorkPointer = term + *term;
4196 else AT.WorkPointer = termout;
4197/*
4198 if ( termflag ) { We have to assign the term to $variable(s)
4199 TermAssign(term);
4200 }
4201*/
4202 AT.NormDepth--;
4203 TermFree(n_llnum,"n_llnum");
4204 TermFree(n_coef,"NormCoef");
4205 return(regval);
4206
4207NormInf:
4208 MLOCK(ErrorMessageLock);
4209 MesPrint("Division by zero during normalization");
4210 MUNLOCK(ErrorMessageLock);
4211 Terminate(-1);
4212
4213NormZZ:
4214 MLOCK(ErrorMessageLock);
4215 MesPrint("0^0 during normalization of term");
4216 MUNLOCK(ErrorMessageLock);
4217 Terminate(-1);
4218
4219NormPRF:
4220 MLOCK(ErrorMessageLock);
4221 MesPrint("0/0 in polyratfun during normalization of term");
4222 MUNLOCK(ErrorMessageLock);
4223 Terminate(-1);
4224
4225NormZero:
4226 *term = 0;
4227 AT.WorkPointer = termout;
4228 AT.NormDepth--;
4229 TermFree(n_llnum,"n_llnum");
4230 TermFree(n_coef,"NormCoef");
4231 return(regval);
4232
4233NormMin:
4234 AT.NormDepth--;
4235 TermFree(n_llnum,"n_llnum");
4236 TermFree(n_coef,"NormCoef");
4237 return(-1);
4238
4239FromNorm:
4240 MLOCK(ErrorMessageLock);
4241 MesCall("Norm");
4242 MUNLOCK(ErrorMessageLock);
4243 AT.NormDepth--;
4244 TermFree(n_llnum,"n_llnum");
4245 TermFree(n_coef,"NormCoef");
4246 return(-1);
4247
4248/*
4249 #] Errors and Finish :
4250*/
4251}
4252
4253/*
4254 #] Normalize :
4255 #[ ExtraSymbol :
4256*/
4257
4258int ExtraSymbol(WORD sym, WORD pow, WORD nsym, WORD *ppsym, WORD *ncoef)
4259{
4260 WORD *m, i;
4261 i = nsym;
4262 m = ppsym - 2;
4263 while ( i > 0 ) {
4264 if ( sym == *m ) {
4265 m++;
4266 if ( pow > 2*MAXPOWER || pow < -2*MAXPOWER
4267 || *m > 2*MAXPOWER || *m < -2*MAXPOWER ) {
4268 MLOCK(ErrorMessageLock);
4269 MesPrint("Illegal wildcard power combination.");
4270 MUNLOCK(ErrorMessageLock);
4271 Terminate(-1);
4272 }
4273 *m += pow;
4274
4275 if ( ( sym <= NumSymbols && sym > -MAXPOWER )
4276 && ( symbols[sym].complex & VARTYPEROOTOFUNITY ) == VARTYPEROOTOFUNITY ) {
4277 *m %= symbols[sym].maxpower;
4278 if ( *m < 0 ) *m += symbols[sym].maxpower;
4279 if ( ( symbols[sym].complex & VARTYPEMINUS ) == VARTYPEMINUS ) {
4280 if ( ( ( symbols[sym].maxpower & 1 ) == 0 ) &&
4281 ( *m >= symbols[sym].maxpower/2 ) ) {
4282 *m -= symbols[sym].maxpower/2; *ncoef = -*ncoef;
4283 }
4284 }
4285 }
4286
4287 if ( *m >= 2*MAXPOWER || *m <= -2*MAXPOWER ) {
4288 MLOCK(ErrorMessageLock);
4289 MesPrint("Power overflow during normalization");
4290 MUNLOCK(ErrorMessageLock);
4291 return(-1);
4292 }
4293 if ( !*m ) {
4294 m--;
4295 while ( i < nsym )
4296 { *m = m[2]; m++; *m = m[2]; m++; i++; }
4297 return(-1);
4298 }
4299 return(0);
4300 }
4301 else if ( sym < *m ) {
4302 m -= 2;
4303 i--;
4304 }
4305 else break;
4306 }
4307 m = ppsym;
4308 while ( i < nsym )
4309 { m--; m[2] = *m; m--; m[2] = *m; i++; }
4310 *m++ = sym;
4311 *m = pow;
4312 return(1);
4313}
4314
4315/*
4316 #] ExtraSymbol :
4317 #[ DoTheta :
4318*/
4319
4320int DoTheta(PHEAD WORD *t)
4321{
4322 GETBIDENTITY
4323 WORD k, *r1, *r2, *tstop, type;
4324 WORD ia, *ta, *tb, *stopa, *stopb;
4325 if ( AC.BracketNormalize ) return(-1);
4326 type = *t;
4327 k = t[1];
4328 tstop = t + k;
4329 t += FUNHEAD;
4330 if ( k <= FUNHEAD ) return(1);
4331 r1 = t;
4332 NEXTARG(r1)
4333 if ( r1 == tstop ) {
4334/*
4335 One argument
4336*/
4337 if ( *t == ARGHEAD ) {
4338 if ( type == THETA ) return(1);
4339 else return(0); /* THETA2 */
4340 }
4341 if ( *t < 0 ) {
4342 if ( *t == -SNUMBER ) {
4343 if ( t[1] < 0 ) return(0);
4344 else {
4345 if ( type == THETA2 && t[1] == 0 ) return(0);
4346 else return(1);
4347 }
4348 }
4349 return(-1);
4350 }
4351 k = t[*t-1];
4352 if ( *t == ABS(k)+1+ARGHEAD ) {
4353 if ( k > 0 ) return(1);
4354 else return(0);
4355 }
4356 return(-1);
4357 }
4358/*
4359 At least two arguments
4360*/
4361 r2 = r1;
4362 NEXTARG(r2)
4363 if ( r2 < tstop ) return(-1); /* More than 2 arguments */
4364/*
4365 Note now that zero has to be treated specially
4366 We take the criteria from the symmetrize routine
4367*/
4368 if ( *t == -SNUMBER && *r1 == -SNUMBER ) {
4369 if ( t[1] > r1[1] ) return(0);
4370 else if ( t[1] < r1[1] ) {
4371 return(1);
4372 }
4373 else if ( type == THETA ) return(1);
4374 else return(0); /* THETA2 */
4375 }
4376 else if ( t[1] == 0 && *t == -SNUMBER ) {
4377 if ( *r1 > 0 ) { }
4378 else if ( *t < *r1 ) return(1);
4379 else if ( *t > *r1 ) return(0);
4380 }
4381 else if ( r1[1] == 0 && *r1 == -SNUMBER ) {
4382 if ( *t > 0 ) { }
4383 else if ( *t < *r1 ) return(1);
4384 else if ( *t > *r1 ) return(0);
4385 }
4386 r2 = AT.WorkPointer;
4387 if ( *t < 0 ) {
4388 ta = r2;
4389 ToGeneral(t,ta,0);
4390 r2 += *r2;
4391 }
4392 else ta = t;
4393 if ( *r1 < 0 ) {
4394 tb = r2;
4395 ToGeneral(r1,tb,0);
4396 }
4397 else tb = r1;
4398 stopa = ta + *ta;
4399 stopb = tb + *tb;
4400 ta += ARGHEAD; tb += ARGHEAD;
4401 while ( ta < stopa ) {
4402 if ( tb >= stopb ) return(0);
4403 if ( ( ia = CompareTerms(BHEAD ta,tb,(WORD)1) ) < 0 ) return(0);
4404 if ( ia > 0 ) return(1);
4405 ta += *ta;
4406 tb += *tb;
4407 }
4408 if ( type == THETA ) return(1);
4409 else return(0); /* THETA2 */
4410}
4411
4412/*
4413 #] DoTheta :
4414 #[ DoDelta :
4415*/
4416
4417int DoDelta(WORD *t)
4418{
4419 WORD k, *r1, *r2, *tstop, isnum, isnum2, type = *t;
4420 if ( AC.BracketNormalize ) return(-1);
4421 k = t[1];
4422 if ( k <= FUNHEAD ) goto argzero;
4423 if ( k == FUNHEAD+ARGHEAD && t[FUNHEAD] == ARGHEAD ) goto argzero;
4424 tstop = t + k;
4425 t += FUNHEAD;
4426 r1 = t;
4427 NEXTARG(r1)
4428 if ( *t < 0 ) {
4429 k = 1;
4430 if ( *t == -SNUMBER ) { isnum = 1; k = t[1]; }
4431 else isnum = 0;
4432 }
4433 else {
4434 k = t[*t-1];
4435 k = ABS(k);
4436 if ( k == *t-ARGHEAD-1 ) isnum = 1;
4437 else isnum = 0;
4438 k = 1;
4439 }
4440 if ( r1 >= tstop ) { /* Single argument */
4441 if ( !isnum ) return(-1);
4442 if ( k == 0 ) goto argzero;
4443 goto argnonzero;
4444 }
4445 r2 = r1;
4446 NEXTARG(r2)
4447 if ( r2 < tstop ) return(-1);
4448 if ( *r1 < 0 ) {
4449 if ( *r1 == -SNUMBER ) { isnum2 = 1; }
4450 else isnum2 = 0;
4451 }
4452 else {
4453 k = r1[*r1-1];
4454 k = ABS(k);
4455 if ( k == *r1-ARGHEAD-1 ) isnum2 = 1;
4456 else isnum2 = 0;
4457 }
4458 if ( isnum != isnum2 ) return(-1);
4459 tstop = r1;
4460 while ( t < tstop && r1 < r2 ) {
4461 if ( *t != *r1 ) {
4462 if ( !isnum ) return(-1);
4463 goto argnonzero;
4464 }
4465 t++; r1++;
4466 }
4467 if ( t != tstop || r1 != r2 ) {
4468 if ( !isnum ) return(-1);
4469 goto argnonzero;
4470 }
4471argzero:
4472 if ( type == DELTA2 ) return(1);
4473 else return(0);
4474argnonzero:
4475 if ( type == DELTA2 ) return(0);
4476 else return(1);
4477}
4478
4479/*
4480 #] DoDelta :
4481 #[ DoRevert :
4482*/
4483
4484void DoRevert(WORD *fun, WORD *tmp)
4485{
4486 WORD *t, *r, *m, *to, *tt, *mm, i, j;
4487 to = fun + fun[1];
4488 r = fun + FUNHEAD;
4489 while ( r < to ) {
4490 if ( *r <= 0 ) {
4491 if ( *r == -REVERSEFUNCTION ) {
4492 m = r; mm = m+1;
4493 while ( mm < to ) *m++ = *mm++;
4494 to--;
4495 (fun[1])--;
4496 fun[2] |= DIRTYSYMFLAG;
4497 }
4498 else if ( *r <= -FUNCTION ) r++;
4499 else {
4500 if ( *r == -INDEX && r[1] < MINSPEC ) *r = -VECTOR;
4501 r += 2;
4502 }
4503 }
4504 else {
4505 if ( ( *r > ARGHEAD )
4506 && ( r[ARGHEAD+1] == REVERSEFUNCTION )
4507 && ( *r == (r[ARGHEAD]+ARGHEAD) )
4508 && ( r[ARGHEAD] == (r[ARGHEAD+2]+4) )
4509 && ( *(r+*r-3) == 1 )
4510 && ( *(r+*r-2) == 1 )
4511 && ( *(r+*r-1) == 3 ) ) {
4512 mm = r;
4513 r += ARGHEAD + 1;
4514 tt = r + r[1];
4515 r += FUNHEAD;
4516 m = tmp;
4517 t = r;
4518 j = 0;
4519 while ( t < tt ) {
4520 NEXTARG(t)
4521 j++;
4522 }
4523 while ( --j >= 0 ) {
4524 i = j;
4525 t = r;
4526 while ( --i >= 0 ) {
4527 NEXTARG(t)
4528 }
4529 if ( *t > 0 ) {
4530 i = *t;
4531 NCOPY(m,t,i);
4532 }
4533 else if ( *t <= -FUNCTION ) *m++ = *t++;
4534 else { *m++ = *t++; *m++ = *t++; }
4535 }
4536 i = WORDDIF(m,tmp);
4537 m = tmp;
4538 t = mm;
4539 r = t + *t;
4540 NCOPY(t,m,i);
4541 m = r;
4542 r = t;
4543 i = WORDDIF(to,m);
4544 NCOPY(t,m,i);
4545 fun[1] = WORDDIF(t,fun);
4546 to = t;
4547 fun[2] |= DIRTYSYMFLAG;
4548 }
4549 else r += *r;
4550 }
4551 }
4552}
4553
4554/*
4555 #] DoRevert :
4556 #] Normalize :
4557 #[ DetCommu :
4558
4559 Determines the number of terms in an expression that contain
4560 noncommuting objects. This can be used to see whether products of
4561 this expression can be evaluated with binomial coefficients.
4562
4563 We don't try to be fancy. If a term contains noncommuting objects
4564 we are not looking whether they can commute with complete other
4565 terms.
4566
4567 If the number gets too large we cut it off.
4568*/
4569
4570#define MAXNUMBEROFNONCOMTERMS 2
4571
4572WORD DetCommu(WORD *terms)
4573{
4574 WORD *t, *tnext, *tstop;
4575 WORD num = 0;
4576 if ( *terms == 0 ) return(0);
4577 if ( terms[*terms] == 0 ) return(0);
4578 t = terms;
4579 while ( *t ) {
4580 tnext = t + *t;
4581 tstop = tnext - ABS(tnext[-1]);
4582 t++;
4583 while ( t < tstop ) {
4584 if ( *t >= FUNCTION ) {
4585 if ( functions[*t-FUNCTION].commute ) {
4586 num++;
4587 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4588 break;
4589 }
4590 }
4591 else if ( *t == SUBEXPRESSION ) {
4592 if ( cbuf[t[4]].CanCommu[t[2]] ) {
4593 num++;
4594 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4595 break;
4596 }
4597 }
4598 else if ( *t == EXPRESSION ) {
4599 num++;
4600 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4601 break;
4602 }
4603 else if ( *t == DOLLAREXPRESSION ) {
4604/*
4605 Technically this is not correct. We have to test first
4606 whether this is DollarLocalCopy (in TFORM) and if so, use the
4607 local version. Anyway, this should be rare to never
4608 occurring because dollars should be replaced.
4609*/
4610 if ( cbuf[AM.dbufnum].CanCommu[t[2]] ) {
4611 num++;
4612 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4613 break;
4614 }
4615 }
4616 t += t[1];
4617 }
4618 t = tnext;
4619 }
4620 return(num);
4621}
4622
4623/*
4624 #] DetCommu :
4625 #[ DoesCommu :
4626
4627 Determines the number of noncommuting objects in a term.
4628 If the number gets too large we cut it off.
4629*/
4630
4631WORD DoesCommu(WORD *term)
4632{
4633 WORD *tstop;
4634 WORD num = 0;
4635 if ( *term == 0 ) return(0);
4636 tstop = term + *term;
4637 tstop = tstop - ABS(tstop[-1]);
4638 term++;
4639 while ( term < tstop ) {
4640 if ( ( *term >= FUNCTION ) && ( functions[*term-FUNCTION].commute ) ) {
4641 num++;
4642 if ( num >= MAXNUMBEROFNONCOMTERMS ) return(num);
4643 }
4644 term += term[1];
4645 }
4646 return(num);
4647}
4648
4649/*
4650 #] DoesCommu :
4651 #[ PolyNormPoly :
4652
4653 Normalizes a polynomial
4654*/
4655
4656#ifdef EVALUATEGCD
4657WORD *PolyNormPoly (PHEAD WORD *Poly) {
4658
4659 GETBIDENTITY;
4660 WORD *buffer = AT.WorkPointer;
4661 WORD *p;
4662 if ( NewSort(BHEAD0) ) { Terminate(-1); }
4663 AR.CompareRoutine = (COMPAREDUMMY)(&CompareSymbols);
4664 while ( *Poly ) {
4665 p = Poly + *Poly;
4666 if ( SymbolNormalize(Poly) < 0 ) return(0);
4667 if ( StoreTerm(BHEAD Poly) ) {
4668 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
4669 LowerSortLevel();
4670 Terminate(-1);
4671 }
4672 Poly = p;
4673 }
4674 if ( EndSort(BHEAD buffer,1) < 0 ) {
4675 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
4676 Terminate(-1);
4677 }
4678 p = buffer;
4679 while ( *p ) p += *p;
4680 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
4681 AT.WorkPointer = p + 1;
4682 return(buffer);
4683}
4684#endif
4685
4686/*
4687 #] PolyNormPoly :
4688 #[ EvaluateGcd :
4689
4690 Try to evaluate the GCDFUNCTION gcd_.
4691 This function can have a number of arguments which can be numbers
4692 and/or polynomials. If there are objects that aren't SYMBOLS or numbers
4693 it cannot work currently.
4694
4695 To make this work properly we have to intervene in proces.c
4696 proces.c: if ( Normalize(BHEAD m) ) {
46971060 proces.c: if ( Normalize(BHEAD r) ) {
46981126?proces.c: if ( Normalize(BHEAD term) ) {
4699 proces.c: if ( Normalize(BHEAD AT.WorkPointer) ) goto PasErr;
47002308!proces.c: if ( ( retnorm = Normalize(BHEAD term) ) != 0 ) {
4701 proces.c: ReNumber(BHEAD term); Normalize(BHEAD term);
4702 proces.c: if ( Normalize(BHEAD v) ) Terminate(-1);
4703 proces.c: if ( Normalize(BHEAD w) ) { LowerSortLevel(); goto PolyCall; }
4704 proces.c: if ( Normalize(BHEAD term) ) goto PolyCall;
4705*/
4706#ifdef EVALUATEGCD
4707
4708WORD *EvaluateGcd(PHEAD WORD *subterm)
4709{
4710 GETBIDENTITY
4711 WORD *oldworkpointer = AT.WorkPointer, *work1, *work2, *work3;
4712 WORD *t, *tt, *ttt, *t1, *t2, *t3, *t4, *tstop;
4713 WORD ct, nnum;
4714 UWORD gcdnum, stor;
4715 WORD *lnum=n_llnum+1;
4716 WORD *num1, *num2, *num3, *den1, *den2, *den3;
4717 WORD sizenum1, sizenum2, sizenum3, sizeden1, sizeden2, sizeden3;
4718 int i, isnumeric = 0, numarg = 0 /*, sizearg */;
4719 LONG size;
4720/*
4721 Step 1: Look for -SNUMBER or -SYMBOL arguments.
4722 If encountered, treat everybody with it.
4723*/
4724 tt = subterm + subterm[1]; t = subterm + FUNHEAD;
4725
4726 while ( t < tt ) {
4727 numarg++;
4728 if ( *t == -SNUMBER ) {
4729 if ( t[1] == 0 ) {
4730gcdzero:;
4731 MLOCK(ErrorMessageLock);
4732 MesPrint("Trying to take the GCD involving a zero term.");
4733 MUNLOCK(ErrorMessageLock);
4734 return(0);
4735 }
4736 gcdnum = ABS(t[1]);
4737 t1 = subterm + FUNHEAD;
4738 while ( gcdnum > 1 && t1 < tt ) {
4739 if ( *t1 == -SNUMBER ) {
4740 stor = ABS(t1[1]);
4741 if ( stor == 0 ) goto gcdzero;
4742 if ( GcdLong(BHEAD (UWORD *)&stor,1,(UWORD *)&gcdnum,1,
4743 (UWORD *)lnum,&nnum) ) goto FromGCD;
4744 gcdnum = lnum[0];
4745 t1 += 2;
4746 continue;
4747 }
4748 else if ( *t1 == -SYMBOL ) goto gcdisone;
4749 else if ( *t1 < 0 ) goto gcdillegal;
4750/*
4751 Now we have to go through all the terms in the argument.
4752 This includes long numbers.
4753*/
4754 ttt = t1 + *t1;
4755 ct = *ttt; *ttt = 0;
4756 if ( t1[1] != 0 ) { /* First normalize the argument */
4757 t1 = PolyNormPoly(BHEAD t1+ARGHEAD);
4758 }
4759 else t1 += ARGHEAD;
4760 while ( *t1 ) {
4761 t1 += *t1;
4762 i = ABS(t1[-1]);
4763 t2 = t1 - i;
4764 i = (i-1)/2;
4765 t3 = t2+i-1;
4766 while ( t3 > t2 && *t3 == 0 ) { t3--; i--; }
4767 if ( GcdLong(BHEAD (UWORD *)t2,(WORD)i,(UWORD *)&gcdnum,1,
4768 (UWORD *)lnum,&nnum) ) {
4769 *ttt = ct;
4770 goto FromGCD;
4771 }
4772 gcdnum = lnum[0];
4773 if ( gcdnum == 1 ) {
4774 *ttt = ct;
4775 goto gcdisone;
4776 }
4777 }
4778 *ttt = ct;
4779 t1 = ttt;
4780 AT.WorkPointer = oldworkpointer;
4781 }
4782 if ( gcdnum == 1 ) goto gcdisone;
4783 oldworkpointer[0] = 4;
4784 oldworkpointer[1] = gcdnum;
4785 oldworkpointer[2] = 1;
4786 oldworkpointer[3] = 3;
4787 oldworkpointer[4] = 0;
4788 AT.WorkPointer = oldworkpointer + 5;
4789 return(oldworkpointer);
4790 }
4791 else if ( *t == -SYMBOL ) {
4792 t1 = subterm + FUNHEAD;
4793 i = t[1];
4794 while ( t1 < tt ) {
4795 if ( *t1 == -SNUMBER ) goto gcdisone;
4796 if ( *t1 == -SYMBOL ) {
4797 if ( t1[1] != i ) goto gcdisone;
4798 t1 += 2; continue;
4799 }
4800 if ( *t1 < 0 ) goto gcdillegal;
4801 ttt = t1 + *t1;
4802 ct = *ttt; *ttt = 0;
4803 if ( t1[1] != 0 ) { /* First normalize the argument */
4804 t2 = PolyNormPoly(BHEAD t1+ARGHEAD);
4805 }
4806 else t2 = t1 + ARGHEAD;
4807 while ( *t2 ) {
4808 t3 = t2+1;
4809 t2 = t2 + *t2;
4810 tstop = t2 - ABS(t2[-1]);
4811 while ( t3 < tstop ) {
4812 if ( *t3 != SYMBOL ) {
4813 *ttt = ct;
4814 goto gcdillegal;
4815 }
4816 t4 = t3 + 2;
4817 t3 += t3[1];
4818 while ( t4 < t3 ) {
4819 if ( *t4 == i && t4[1] > 0 ) goto nextterminarg;
4820 t4 += 2;
4821 }
4822 }
4823 *ttt = ct;
4824 goto gcdisone;
4825nextterminarg:;
4826 }
4827 *ttt = ct;
4828 t1 = ttt;
4829 AT.WorkPointer = oldworkpointer;
4830 }
4831 oldworkpointer[0] = 8;
4832 oldworkpointer[1] = SYMBOL;
4833 oldworkpointer[2] = 4;
4834 oldworkpointer[3] = t[1];
4835 oldworkpointer[4] = 1;
4836 oldworkpointer[5] = 1;
4837 oldworkpointer[6] = 1;
4838 oldworkpointer[7] = 3;
4839 oldworkpointer[8] = 0;
4840 AT.WorkPointer = oldworkpointer+9;
4841 return(oldworkpointer);
4842 }
4843 else if ( *t < 0 ) {
4844gcdillegal:;
4845 MLOCK(ErrorMessageLock);
4846 MesPrint("Illegal object in gcd_ function. Object not a number or a symbol.");
4847 MUNLOCK(ErrorMessageLock);
4848 goto FromGCD;
4849 }
4850 else if ( ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4851 else if ( t[1] != 0 ) {
4852 ttt = t + *t; ct = *ttt; *ttt = 0;
4853 t = PolyNormPoly(BHEAD t+ARGHEAD);
4854 *ttt = ct;
4855 if ( t[*t] == 0 && ABS(t[*t-1]) == *t-ARGHEAD-1 ) isnumeric = numarg;
4856 AT.WorkPointer = oldworkpointer;
4857 t = ttt;
4858 }
4859 t += *t;
4860 }
4861/*
4862 At this point there are only generic arguments.
4863 There are however still two cases:
4864 1: There is an argument that is purely numerical
4865 In that case we have to take the gcd of all coefficients
4866 2: All arguments are nontrivial polynomials.
4867 Here we don't worry so much about the factor. (???)
4868 We know whether case 1 occurs when isnumeric > 0.
4869 We can look up numarg to get a good starting value.
4870*/
4871 AT.WorkPointer = oldworkpointer;
4872 if ( isnumeric ) {
4873 t = subterm + FUNHEAD;
4874 for ( i = 1; i < isnumeric; i++ ) {
4875 NEXTARG(t);
4876 }
4877 if ( t[1] != 0 ) { /* First normalize the argument */
4878 ttt = t + *t; ct = *ttt; *ttt = 0;
4879 t = PolyNormPoly(BHEAD t+ARGHEAD);
4880 *ttt = ct;
4881 }
4882 t += *t;
4883 i = (ABS(t[-1])-1)/2;
4884 den1 = t - 1 - i;
4885 num1 = den1 - i;
4886 sizenum1 = sizeden1 = i;
4887 while ( sizenum1 > 1 && num1[sizenum1-1] == 0 ) sizenum1--;
4888 while ( sizeden1 > 1 && den1[sizeden1-1] == 0 ) sizeden1--;
4889 work1 = AT.WorkPointer+1; work2 = work1+sizenum1;
4890 for ( i = 0; i < sizenum1; i++ ) work1[i] = num1[i];
4891 for ( i = 0; i < sizeden1; i++ ) work2[i] = den1[i];
4892 num1 = work1; den1 = work2;
4893 AT.WorkPointer = work2 = work2 + sizeden1;
4894 t = subterm + FUNHEAD;
4895 while ( t < tt ) {
4896 ttt = t + *t; ct = *ttt; *ttt = 0;
4897 if ( t[1] != 0 ) {
4898 t = PolyNormPoly(BHEAD t+ARGHEAD);
4899 }
4900 else t += ARGHEAD;
4901 while ( *t ) {
4902 t += *t;
4903 i = (ABS(t[-1])-1)/2;
4904 den2 = t - 1 - i;
4905 num2 = den2 - i;
4906 sizenum2 = sizeden2 = i;
4907 while ( sizenum2 > 1 && num2[sizenum2-1] == 0 ) sizenum2--;
4908 while ( sizeden2 > 1 && den2[sizeden2-1] == 0 ) sizeden2--;
4909 num3 = AT.WorkPointer;
4910 if ( GcdLong(BHEAD (UWORD *)num2,sizenum2,(UWORD *)num1,sizenum1,
4911 (UWORD *)num3,&sizenum3) ) goto FromGCD;
4912 sizenum1 = sizenum3;
4913 for ( i = 0; i < sizenum1; i++ ) num1[i] = num3[i];
4914 den3 = AT.WorkPointer;
4915 if ( GcdLong(BHEAD (UWORD *)den2,sizeden2,(UWORD *)den1,sizeden1,
4916 (UWORD *)den3,&sizeden3) ) goto FromGCD;
4917 sizeden1 = sizeden3;
4918 for ( i = 0; i < sizeden1; i++ ) den1[i] = den3[i];
4919 if ( sizenum1 == 1 && num1[0] == 1 && sizeden1 == 1 && den1[1] == 1 )
4920 goto gcdisone;
4921 }
4922 *ttt = ct;
4923 t = ttt;
4924 AT.WorkPointer = work2;
4925 }
4926 AT.WorkPointer = oldworkpointer;
4927/*
4928 Now copy the GCD to the 'output'
4929*/
4930 if ( sizenum1 > sizeden1 ) {
4931 while ( sizenum1 > sizeden1 ) den1[sizeden1++] = 0;
4932 }
4933 else if ( sizenum1 < sizeden1 ) {
4934 while ( sizenum1 < sizeden1 ) num1[sizenum1++] = 0;
4935 }
4936 t = oldworkpointer;
4937 i = 2*sizenum1+1;
4938 *t++ = i+1;
4939 if ( num1 != t ) { NCOPY(t,num1,sizenum1); }
4940 else t += sizenum1;
4941 if ( den1 != t ) { NCOPY(t,den1,sizeden1); }
4942 else t += sizeden1;
4943 *t++ = i;
4944 *t++ = 0;
4945 AT.WorkPointer = t;
4946 return(oldworkpointer);
4947 }
4948/*
4949 Now the real stuff with only polynomials.
4950 Pick up the shortest term to start.
4951 We are a bit brutish about this.
4952*/
4953 t = subterm + FUNHEAD;
4954 AT.WorkPointer += AM.MaxTer/sizeof(WORD);
4955 work2 = AT.WorkPointer;
4956/*
4957 sizearg = subterm[1];
4958*/
4959 i = 0; work3 = 0;
4960 while ( t < tt ) {
4961 i++;
4962 work1 = AT.WorkPointer;
4963 ttt = t + *t; ct = *ttt; *ttt = 0;
4964 t = PolyNormPoly(BHEAD t+ARGHEAD);
4965 if ( *work1 < AT.WorkPointer-work1 ) {
4966/*
4967 sizearg = AT.WorkPointer-work1;
4968*/
4969 numarg = i;
4970 work3 = work1;
4971 }
4972 *ttt = ct; t = ttt;
4973 }
4974 *AT.WorkPointer++ = 0;
4975/*
4976 We have properly normalized arguments and the shortest is indicated in work3
4977*/
4978 work1 = work3;
4979 while ( *work2 ) {
4980 if ( work2 != work3 ) {
4981 work1 = PolyGCD2(BHEAD work1,work2);
4982 }
4983 while ( *work2 ) work2 += *work2;
4984 work2++;
4985 }
4986 work2 = work1;
4987 while ( *work2 ) work2 += *work2;
4988 size = work2 - work1 + 1;
4989 t = oldworkpointer;
4990 NCOPY(t,work1,size);
4991 AT.WorkPointer = t;
4992 return(oldworkpointer);
4993
4994gcdisone:;
4995 oldworkpointer[0] = 4;
4996 oldworkpointer[1] = 1;
4997 oldworkpointer[2] = 1;
4998 oldworkpointer[3] = 3;
4999 oldworkpointer[4] = 0;
5000 AT.WorkPointer = oldworkpointer+5;
5001 return(oldworkpointer);
5002FromGCD:
5003 MLOCK(ErrorMessageLock);
5004 MesCall("EvaluateGcd");
5005 MUNLOCK(ErrorMessageLock);
5006 return(0);
5007}
5008
5009#endif
5010
5011/*
5012 #] EvaluateGcd :
5013 #[ TreatPolyRatFun :
5014
5015 if ( AR.PolyFunExp == 1 ) we have to trim the contents of the polyratfun
5016 down to its most divergent term and give it coefficient +1. This is done
5017 by taking the terms with the least power in the variable in the numerator
5018 and in the denominator and then combine them.
5019 Answer is either PolyRatFun(ep^n,1) or PolyRatFun(1,1) or PolyRatFun(1,ep^n)
5020*/
5021
5022int TreatPolyRatFun(PHEAD WORD *prf)
5023{
5024 WORD *t, *tstop, *r, *rstop, *m, *mstop;
5025 WORD exp1 = MAXPOWER, exp2 = MAXPOWER;
5026 t = prf+FUNHEAD;
5027 if ( *t < 0 ) {
5028 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
5029 if ( exp1 > 1 ) exp1 = 1;
5030 t += 2;
5031 }
5032 else {
5033 if ( exp1 > 0 ) exp1 = 0;
5034 NEXTARG(t)
5035 }
5036 }
5037 else {
5038 tstop = t + *t;
5039 t += ARGHEAD;
5040 while ( t < tstop ) {
5041/*
5042 Now look for the minimum power of AR.PolyFunVar
5043*/
5044 r = t+1;
5045 t += *t;
5046 rstop = t - ABS(t[-1]);
5047 while ( r < rstop ) {
5048 if ( *r != SYMBOL ) { r += r[1]; continue; }
5049 m = r;
5050 mstop = m + m[1];
5051 m += 2;
5052 while ( m < mstop ) {
5053 if ( *m == AR.PolyFunVar ) {
5054 if ( m[1] < exp1 ) exp1 = m[1];
5055 break;
5056 }
5057 m += 2;
5058 }
5059 if ( m == mstop ) {
5060 if ( exp1 > 0 ) exp1 = 0;
5061 }
5062 break;
5063 }
5064 if ( r == rstop ) {
5065 if ( exp1 > 0 ) exp1 = 0;
5066 }
5067 }
5068 t = tstop;
5069 }
5070 if ( *t < 0 ) {
5071 if ( *t == -SYMBOL && t[1] == AR.PolyFunVar ) {
5072 if ( exp2 > 1 ) exp2 = 1;
5073 }
5074 else {
5075 if ( exp2 > 0 ) exp2 = 0;
5076 }
5077 }
5078 else {
5079 tstop = t + *t;
5080 t += ARGHEAD;
5081 while ( t < tstop ) {
5082/*
5083 Now look for the minimum power of AR.PolyFunVar
5084*/
5085 r = t+1;
5086 t += *t;
5087 rstop = t - ABS(t[-1]);
5088 while ( r < rstop ) {
5089 if ( *r != SYMBOL ) { r += r[1]; continue; }
5090 m = r;
5091 mstop = m + m[1];
5092 m += 2;
5093 while ( m < mstop ) {
5094 if ( *m == AR.PolyFunVar ) {
5095 if ( m[1] < exp2 ) exp2 = m[1];
5096 break;
5097 }
5098 m += 2;
5099 }
5100 if ( m == mstop ) {
5101 if ( exp2 > 0 ) exp2 = 0;
5102 }
5103 break;
5104 }
5105 if ( r == rstop ) {
5106 if ( exp2 > 0 ) exp2 = 0;
5107 }
5108 }
5109 }
5110/*
5111 Now we can compose the output.
5112 Notice that the output can never be longer than the input provided
5113 we never can have arguments that consist of just a function.
5114*/
5115 exp1 = exp1-exp2;
5116/* if ( exp1 > 0 ) exp1 = 0; */
5117 t = prf+FUNHEAD;
5118 if ( exp1 == 0 ) {
5119 *t++ = -SNUMBER; *t++ = 1;
5120 *t++ = -SNUMBER; *t++ = 1;
5121 }
5122 else if ( exp1 > 0 ) {
5123 if ( exp1 == 1 ) {
5124 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
5125 }
5126 else {
5127 *t++ = 8+ARGHEAD;
5128 *t++ = 0;
5129 FILLARG(t);
5130 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
5131 *t++ = exp1; *t++ = 1; *t++ = 1; *t++ = 3;
5132 }
5133 *t++ = -SNUMBER; *t++ = 1;
5134 }
5135 else {
5136 *t++ = -SNUMBER; *t++ = 1;
5137 if ( exp1 == -1 ) {
5138 *t++ = -SYMBOL; *t++ = AR.PolyFunVar;
5139 }
5140 else {
5141 *t++ = 8+ARGHEAD;
5142 *t++ = 0;
5143 FILLARG(t);
5144 *t++ = 8; *t++ = SYMBOL; *t++ = 4; *t++ = AR.PolyFunVar;
5145 *t++ = -exp1; *t++ = 1; *t++ = 1; *t++ = 3;
5146 }
5147 }
5148 prf[2] = 0; /* Clean */
5149 prf[1] = t - prf;
5150 return(0);
5151}
5152
5153/*
5154 #] TreatPolyRatFun :
5155 #[ DropCoefficient :
5156*/
5157
5158void DropCoefficient(PHEAD WORD *term)
5159{
5160 GETBIDENTITY
5161 WORD *t = term + *term;
5162 WORD n, na;
5163 n = t[-1]; na = ABS(n);
5164 t -= na;
5165 if ( n == 3 && t[0] == 1 && t[1] == 1 ) return;
5166 *AN.RepPoint = 1;
5167 t[0] = 1; t[1] = 1; t[2] = 3;
5168 *term -= (na-3);
5169}
5170
5171/*
5172 #] DropCoefficient :
5173 #[ DropSymbols :
5174*/
5175
5176void DropSymbols(PHEAD WORD *term)
5177{
5178 GETBIDENTITY
5179 WORD *tend = term + *term, *t1, *t2, *tstop;
5180 tstop = tend - ABS(tend[-1]);
5181 t1 = term+1;
5182 while ( t1 < tstop ) {
5183 if ( *t1 == SYMBOL ) {
5184 *AN.RepPoint = 1;
5185 t2 = t1+t1[1];
5186 while ( t2 < tend ) *t1++ = *t2++;
5187 *term = t1 - term;
5188 break;
5189 }
5190 t1 += t1[1];
5191 }
5192}
5193
5194/*
5195 #] DropSymbols :
5196 #[ SymbolNormalize :
5197*/
5206int SymbolNormalize(WORD *term)
5207{
5208 GETIDENTITY
5209 WORD *t, *b, *bb, *tt, *m, *tstop;
5210 // Here we use a stack-allocated array, since things are much smaller
5211 // compared to the full Normalize routine.
5212 WORD buffer[7*NORMSIZE];
5213 int i;
5214 b = buffer;
5215 *b++ = SYMBOL; *b++ = 2;
5216 t = term + *term;
5217 tstop = t - ABS(t[-1]);
5218 t = term + 1;
5219 while ( t < tstop ) { /* Step 1: collect symbols */
5220 if ( *t == SYMBOL && t < tstop ) {
5221 for ( i = 2; i < t[1]; i += 2 ) {
5222 const WORD sym = t[i];
5223 const WORD pow = t[i+1];
5224 bb = buffer+2;
5225 while ( bb < b ) {
5226 if ( bb[0] == sym ) { /* add powers */
5227 bb[1] += pow;
5228 if ( bb[1] > MAXPOWER || bb[1] < -MAXPOWER ) {
5229 MLOCK(ErrorMessageLock);
5230 MesPrint("Power in SymbolNormalize out of range");
5231 MUNLOCK(ErrorMessageLock);
5232 return(-1);
5233 }
5234 if ( bb[1] == 0 ) {
5235 b -= 2;
5236 while ( bb < b ) {
5237 bb[0] = bb[2]; bb[1] = bb[3]; bb += 2;
5238 }
5239 }
5240 goto Nexti;
5241 }
5242 else if ( bb[0] > sym ) { /* insert it */
5243 m = b;
5244 while ( m > bb ) { m[1] = m[-1]; m[0] = m[-2]; m -= 2; }
5245 b += 2;
5246 bb[0] = sym;
5247 bb[1] = pow;
5248 goto Nexti;
5249 }
5250 bb += 2;
5251 }
5252 if ( bb >= b ) { /* add it to the end */
5253 *b++ = sym; *b++ = pow;
5254 }
5255Nexti:;
5256 }
5257 }
5258 else {
5259 MLOCK(ErrorMessageLock);
5260 MesPrint("Illegal term in SymbolNormalize");
5261 MUNLOCK(ErrorMessageLock);
5262 return(-1);
5263 }
5264 t += t[1];
5265 }
5266 buffer[1] = b - buffer;
5267/*
5268 Veto negative powers
5269*/
5270 if ( AT.LeaveNegative == 0 ) {
5271 b = buffer; bb = b + b[1]; b += 3;
5272 while ( b < bb ) {
5273 if ( *b < 0 ) {
5274 MLOCK(ErrorMessageLock);
5275 MesPrint("Negative power in SymbolNormalize");
5276 MUNLOCK(ErrorMessageLock);
5277 return(-1);
5278 }
5279 b += 2;
5280 }
5281 }
5282/*
5283 Now we use the fact that the new term will not be longer than the old one
5284 Actually it should be shorter when there is more than one subterm!
5285 Copy back.
5286*/
5287 i = buffer[1];
5288 b = buffer; tt = term + 1;
5289 if ( i > 2 ) { NCOPY(tt,b,i) }
5290 if ( tt < tstop ) {
5291 i = term[*term-1];
5292 if ( i < 0 ) i = -i;
5293 *term -= (tstop-tt);
5294 NCOPY(tt,tstop,i)
5295 }
5296 return(0);
5297}
5298
5299/*
5300 #] SymbolNormalize :
5301 #[ TestFunFlag :
5302
5303 Tests whether a function still has unsubstituted subexpressions
5304 This function has its dirtyflag on!
5305*/
5306
5307int TestFunFlag(PHEAD WORD *tfun)
5308{
5309 WORD *t, *tstop, *r, *rstop, *m, *mstop;
5310 if ( functions[*tfun-FUNCTION].spec <= 0 ) return(0);
5311 tstop = tfun + tfun[1];
5312 t = tfun + FUNHEAD;
5313 while ( t < tstop ) {
5314 if ( *t < 0 ) { NEXTARG(t); continue; }
5315 rstop = t + *t;
5316 if ( t[1] == 0 ) { t = rstop; continue; }
5317 r = t + ARGHEAD;
5318 while ( r < rstop ) { /* Here we loop over terms */
5319 m = r+1; mstop = r+*r; mstop -= ABS(mstop[-1]);
5320 while ( m < mstop ) { /* Loop over the subterms */
5321 if ( *m == SUBEXPRESSION || *m == EXPRESSION || *m == DOLLAREXPRESSION ) return(1);
5322 if ( ( *m >= FUNCTION ) && ( ( m[2] & DIRTYFLAG ) == DIRTYFLAG )
5323 && ( *m != REPLACEMENT ) && TestFunFlag(BHEAD m) ) return(1);
5324 m += m[1];
5325 }
5326 r += *r;
5327 }
5328 t += *t;
5329 }
5330 return(0);
5331}
5332
5333/*
5334 #] TestFunFlag :
5335 #[ BracketNormalize :
5336*/
5337
5338int BracketNormalize(PHEAD WORD *term)
5339{
5340 WORD *stop = term+*term-3, *t, *tt, *tstart, *r;
5341 WORD *oldwork = AT.WorkPointer;
5342 WORD *termout;
5343 WORD i, ii, j;
5344 termout = AT.WorkPointer = term+*term;
5345/*
5346 First collect all functions and sort them
5347*/
5348 tt = termout+1; t = term+1;
5349 while ( t < stop ) {
5350 if ( *t >= FUNCTION ) { i = t[1]; NCOPY(tt,t,i); }
5351 else t += t[1];
5352 }
5353 if ( tt > termout+1 && tt-termout-1 > termout[2] ) { /* sorting */
5354 r = termout+1; ii = tt-r;
5355 for ( i = 0; i < ii-FUNHEAD; i += FUNHEAD ) { /* Bubble sort */
5356 for ( j = i+FUNHEAD; j > 0; j -= FUNHEAD ) {
5357 if ( functions[r[j-FUNHEAD]-FUNCTION].commute
5358 && functions[r[j]-FUNCTION].commute == 0 ) break;
5359 if ( r[j-FUNHEAD] > r[j] ) EXCH(r[j-FUNHEAD],r[j])
5360 else break;
5361 }
5362 }
5363 }
5364
5365 tstart = tt; t = term + 1; *tt++ = DELTA; *tt++ = 2;
5366 while ( t < stop ) {
5367 if ( *t == DELTA ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5368 else t += t[1];
5369 }
5370 if ( tstart[1] > 2 ) {
5371 for ( r = tstart+2; r < tstart+tstart[1]; r += 2 ) {
5372 if ( r[0] > r[1] ) EXCH(r[0],r[1])
5373 }
5374 }
5375 if ( tstart[1] > 4 ) { /* sorting */
5376 r = tstart+2; ii = tstart[1]-2;
5377 for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5378 for ( j = i+2; j > 0; j -= 2 ) {
5379 if ( r[j-2] > r[j] ) {
5380 EXCH(r[j-2],r[j])
5381 EXCH(r[j-1],r[j+1])
5382 }
5383 else if ( r[j-2] < r[j] ) break;
5384 else {
5385 if ( r[j-1] > r[j+1] ) EXCH(r[j-1],r[j+1])
5386 else break;
5387 }
5388 }
5389 }
5390 tt = tstart+tstart[1];
5391 }
5392 else if ( tstart[1] == 2 ) { tt = tstart; }
5393 else tt = tstart+4;
5394
5395 tstart = tt; t = term + 1; *tt++ = INDEX; *tt++ = 2;
5396 while ( t < stop ) {
5397 if ( *t == INDEX ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5398 else t += t[1];
5399 }
5400 if ( tstart[1] >= 4 ) { /* sorting */
5401 r = tstart+2; ii = tstart[1]-2;
5402 for ( i = 0; i < ii-1; i += 1 ) { /* Bubble sort */
5403 for ( j = i+1; j > 0; j -= 1 ) {
5404 if ( r[j-1] > r[j] ) EXCH(r[j-1],r[j])
5405 else break;
5406 }
5407 }
5408 tt = tstart+tstart[1];
5409 }
5410 else if ( tstart[1] == 2 ) { tt = tstart; }
5411 else tt = tstart+3;
5412
5413 tstart = tt; t = term + 1; *tt++ = DOTPRODUCT; *tt++ = 2;
5414 while ( t < stop ) {
5415 if ( *t == DOTPRODUCT ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5416 else t += t[1];
5417 }
5418 if ( tstart[1] > 5 ) { /* sorting */
5419 r = tstart+2; ii = tstart[1]-2;
5420 for ( i = 0; i < ii; i += 3 ) {
5421 if ( r[i] > r[i+1] ) EXCH(r[i],r[i+1])
5422 }
5423 for ( i = 0; i < ii-3; i += 3 ) { /* Bubble sort */
5424 for ( j = i+3; j > 0; j -= 3 ) {
5425 if ( r[j-3] < r[j] ) break;
5426 if ( r[j-3] > r[j] ) {
5427 EXCH(r[j-3],r[j])
5428 EXCH(r[j-2],r[j+1])
5429 }
5430 else {
5431 if ( r[j-2] > r[j+1] ) EXCH(r[j-2],r[j+1])
5432 else break;
5433 }
5434 }
5435 }
5436 tt = tstart+tstart[1];
5437 }
5438 else if ( tstart[1] == 2 ) { tt = tstart; }
5439 else {
5440 if ( tstart[2] > tstart[3] ) EXCH(tstart[2],tstart[3])
5441 tt = tstart+5;
5442 }
5443
5444 tstart = tt; t = term + 1; *tt++ = SYMBOL; *tt++ = 2;
5445 while ( t < stop ) {
5446 if ( *t == SYMBOL ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5447 else t += t[1];
5448 }
5449 if ( tstart[1] > 4 ) { /* sorting */
5450 r = tstart+2; ii = tstart[1]-2;
5451 for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5452 for ( j = i+2; j > 0; j -= 2 ) {
5453 if ( r[j-2] > r[j] ) EXCH(r[j-2],r[j])
5454 else break;
5455 }
5456 }
5457 tt = tstart+tstart[1];
5458 }
5459 else if ( tstart[1] == 2 ) { tt = tstart; }
5460 else tt = tstart+4;
5461
5462 tstart = tt; t = term + 1; *tt++ = SETSET; *tt++ = 2;
5463 while ( t < stop ) {
5464 if ( *t == SETSET ) { i = t[1]-2; t += 2; tstart[1] += i; NCOPY(tt,t,i); }
5465 else t += t[1];
5466 }
5467 if ( tstart[1] > 4 ) { /* sorting */
5468 r = tstart+2; ii = tstart[1]-2;
5469 for ( i = 0; i < ii-2; i += 2 ) { /* Bubble sort */
5470 for ( j = i+2; j > 0; j -= 2 ) {
5471 if ( r[j-2] > r[j] ) {
5472 EXCH(r[j-2],r[j])
5473 EXCH(r[j-1],r[j+1])
5474 }
5475 else break;
5476 }
5477 }
5478 tt = tstart+tstart[1];
5479 }
5480 else if ( tstart[1] == 2 ) { tt = tstart; }
5481 else tt = tstart+4;
5482 *tt++ = 1; *tt++ = 1; *tt++ = 3;
5483 t = term; i = *termout = tt - termout; tt = termout;
5484 NCOPY(t,tt,i);
5485 AT.WorkPointer = oldwork;
5486 return(0);
5487}
5488
5489/*
5490 #] BracketNormalize :
5491*/
AddRHS
WORD * AddRHS(int num, int type)
Definition comtool.c:214
DoubleCbuffer
WORD * DoubleCbuffer(int num, WORD *w, int par)
Definition comtool.c:143
GetFirstTerm
int GetFirstTerm(WORD *, int, int)
Definition execute.c:1966
CompCoef
WORD CompCoef(WORD *, WORD *)
Definition reken.c:3048
EndSort
LONG EndSort(PHEAD WORD *, int)
Definition sort.c:486
LowerSortLevel
void LowerSortLevel(void)
Definition sort.c:4703
StoreTerm
int StoreTerm(PHEAD WORD *)
Definition sort.c:4285
NextPrime
WORD NextPrime(PHEAD WORD)
Definition reken.c:3665
NewSort
int NewSort(PHEAD0)
Definition sort.c:395
Compare1
WORD Compare1(PHEAD WORD *, WORD *, WORD)
Definition sort.c:2375
CompareSymbols
WORD CompareSymbols(PHEAD WORD *, WORD *, WORD)
Definition sort.c:2838
SymbolNormalize
int SymbolNormalize(WORD *term)
Definition normal.c:5206
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