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