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