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