1 /** @file notation.c
2 *
3 * Contains the functions that deal with the rewriting and manipulation
4 * of expressions/terms in polynomial representation.
5 */
6 /* #[ License : */
7 /*
8 * Copyright (C) 1984-2017 J.A.M. Vermaseren
9 * When using this file you are requested to refer to the publication
10 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
11 * This is considered a matter of courtesy as the development was paid
12 * for by FOM the Dutch physics granting agency and we would like to
13 * be able to track its scientific use to convince FOM of its value
14 * for the community.
15 *
16 * This file is part of FORM.
17 *
18 * FORM is free software: you can redistribute it and/or modify it under the
19 * terms of the GNU General Public License as published by the Free Software
20 * Foundation, either version 3 of the License, or (at your option) any later
21 * version.
22 *
23 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
24 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
25 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
26 * details.
27 *
28 * You should have received a copy of the GNU General Public License along
29 * with FORM. If not, see <http://www.gnu.org/licenses/>.
30 */
31 /* #] License : */
32 /*
33 #[ Includes :
34 */
35
36 #include "form3.h"
37
38 /*
39 #] Includes :
40 #[ NormPolyTerm :
41
42 Brings a term to normal form.
43
44 This routine knows objects of the following types:
45 SYMBOL
46 HAAKJE
47 SNUMBER
48 LNUMBER
49 The SNUMBER and LNUMBER are worked into the coefficient.
50 One of the essences here is that everything can be done in place.
51 */
52
NormPolyTerm(PHEAD WORD * term)53 int NormPolyTerm(PHEAD WORD *term)
54 {
55 WORD *tcoef, ncoef, *tstop, *tfill, *t, *tt;
56 int equal, i;
57 WORD *r1, *r2, *r3, *r4, *r5, *rfirst, rv;
58 WORD *lnum, nnum; /* Scratch, originally for factorials */
59 /*
60 One: find the coefficient
61 */
62 tcoef = term+*term;
63 ncoef = tcoef[-1];
64 tstop = tcoef - ABS(tcoef[-1]);
65 tfill = t = term + 1;
66 rfirst = 0;
67 if ( t >= tstop ) { return(*term); }
68 while ( t < tstop ) {
69 switch ( *t ) {
70 case SYMBOL:
71 if ( rfirst == 0 ) {
72 /*
73 Here we only need to sort
74 1: assume no equals. Bubble.
75 */
76 rfirst = t;
77 r2 = rfirst+4; tt = r3 = t + t[1]; equal = 0;
78 while ( r2 < r3 ) {
79 r1 = r2 - 2;
80 if ( *r2 > *r1 ) { r2 += 2; continue; }
81 if ( *r2 == *r1 ) { r2 += 2; equal = 1; continue; }
82 rv = *r1; *r1 = *r2; *r2 = rv;
83 r1 -= 2; r2 -= 2; r4 = r2 + 2;
84 while ( r1 > t ) {
85 if ( *r2 >= *r1 ) { r2 = r4; break; }
86 rv = *r1; *r1 = *r2; *r2 = rv;
87 r1 -= 2; r2 -= 2;
88 }
89 }
90 /*
91 2: hunt down the equal objects
92 postpone eliminating zero powers.
93 */
94 if ( equal ) {
95 r1 = t+2; r2 = r1+2;
96 while ( r2 < r3 ) {
97 if ( *r1 == *r2 ) {
98 r1[1] += r2[1];
99 r4 = r2+2;
100 while ( r4 < r3 ) *r2++ = *r4++;
101 t[1] -= 2;
102 r2 = r1 + 2; r3 -= 2;
103 }
104 }
105 }
106 }
107 else {
108 /*
109 Here we only need to insert
110 */
111 r1 = t + 2; tt = r3 = t + t[1];
112 while ( r1 < r3 ) {
113 r2 = rfirst+2; r4 = rfirst + rfirst[1];
114 while ( r2 < r4 ) {
115 if ( *r1 == *r2 ) {
116 r2[1] += r1[1];
117 break;
118 }
119 else if ( *r2 > *r1 ) {
120 r5 = r4;
121 while ( r5 > r2 ) { r5[1] = r5[-1]; r5[0] = r5[-2]; r5 -= 2; }
122 rfirst[1] += 2;
123 *r2 = *r1; r2[1] = r1[1];
124 break;
125 }
126 r2 += 2;
127 }
128 if ( r2 == r4 ) {
129 rfirst[1] += 2;
130 *r2++ = *r1++; *r2++ = *r1++;
131 }
132 else r1 += 2;
133 }
134 }
135 t = tt;
136 break;
137 case HAAKJE: /* Here we skip brackets */
138 t += t[1];
139 break;
140 case SNUMBER:
141 if ( t[2] < 0 ) {
142 t[2] = -t[2];
143 if ( t[3] & 1 ) ncoef = -ncoef;
144 }
145 else if ( t[2] == 0 ) {
146 if ( t[3] < 0 ) goto NormInf;
147 goto NormZero;
148 }
149 lnum = TermMalloc("lnum");
150 lnum[0] = t[2];
151 nnum = 1;
152 if ( t[3] && RaisPow(BHEAD (UWORD *)lnum,&nnum,(UWORD)(ABS(t[3]))) ) goto FromNorm;
153 ncoef = REDLENG(ncoef);
154 if ( t[3] < 0 ) {
155 if ( Divvy(BHEAD (UWORD *)tstop,&ncoef,(UWORD *)lnum,nnum) )
156 goto FromNorm;
157 }
158 else if ( t[3] > 0 ) {
159 if ( Mully(BHEAD (UWORD *)tstop,&ncoef,(UWORD *)lnum,nnum) )
160 goto FromNorm;
161 }
162 ncoef = INCLENG(ncoef);
163 t += t[1];
164 TermFree(lnum,"lnum");
165 break;
166 case LNUMBER:
167 ncoef = REDLENG(ncoef);
168 if ( Mully(BHEAD (UWORD *)tstop,&ncoef,(UWORD *)(t+3),t[2]) ) goto FromNorm;
169 ncoef = INCLENG(ncoef);
170 t += t[1];
171 break;
172 default:
173 MLOCK(ErrorMessageLock);
174 MesPrint("Illegal code in NormPolyTerm");
175 MUNLOCK(ErrorMessageLock);
176 Terminate(-1);
177 break;
178 }
179 }
180 /*
181 Now we try to eliminate objects to the power zero.
182 */
183 if ( rfirst ) {
184 r2 = rfirst+2;
185 r3 = rfirst + rfirst[1];
186 while ( r2 < r3 ) {
187 if ( r2[1] == 0 ) {
188 r1 = r2 + 2;
189 while ( r1 < r3 ) { r1[-2] = r1[0]; r1[-1] = r1[1]; r1 += 2; }
190 r3 -= 2;
191 rfirst[1] -= 2;
192 }
193 else { r2 += 2; }
194 }
195 if ( rfirst[1] < 4 ) rfirst = 0;
196 }
197 /*
198 Finally we put the term together
199 */
200 if ( rfirst ) {
201 i = rfirst[1];
202 NCOPY(tfill,rfirst,i)
203 }
204 i = ABS(ncoef)-1;
205 NCOPY(tfill,tstop,i)
206 *tfill++ = ncoef;
207 *term = tfill - term;
208 return(*term);
209 NormZero:
210 *term = 0;
211 return(0);
212 NormInf:
213 MLOCK(ErrorMessageLock);
214 MesPrint("0^0 in NormPolyTerm");
215 MUNLOCK(ErrorMessageLock);
216 Terminate(-1);
217 return(-1);
218 FromNorm:
219 MLOCK(ErrorMessageLock);
220 MesCall("NormPolyTerm");
221 MUNLOCK(ErrorMessageLock);
222 Terminate(-1);
223 return(-1);
224 }
225
226 /*
227 #] NormPolyTerm :
228 #[ ComparePoly :
229 */
230 /**
231 * Compares two terms. The answer is:
232 * 0 equal ( with exception of the coefficient )
233 * >0 term1 comes first.
234 * <0 term2 comes first.
235 * This routine should not return an error condition.
236 *
237 * The address of this routine is to be put in AR.CompareRoutine when we
238 * want to use it for sorting.
239 * This makes all existing code work properly and we can just replace the
240 * routine on a thread by thread basis (each thread has its own AR struct).
241 * Don't forget to put the old routine back afterwards!
242 *
243 * @param term1 First input term
244 * @param term2 Second input term
245 * @param level Not used for polynomials
246 * @return 0 equal ( with exception of the coefficient if level == 0. )
247 * >0 term1 comes first.
248 * <0 term2 comes first.
249 */
250
251 #ifdef WITHCOMPAREPOLY
252
ComparePoly(WORD * term1,WORD * term2,WORD level)253 WORD ComparePoly(WORD *term1, WORD *term2, WORD level)
254 {
255 WORD *t1, *t2, *t3, *t4, *tstop1, *tstop2;
256 tstop1 = term1 + *term1;
257 tstop1 -= ABS(tstop1[-1]);
258 tstop2 = term2 + *term2;
259 tstop2 -= ABS(tstop2[-1]);
260 t1 = term1+1;
261 t2 = term2+1;
262 while ( t1 < tstop1 && t2 < tstop2 ) {
263 if ( *t1 == *t2 ) {
264 if ( *t1 == HAAKJE ) {
265 if ( t1[2] != t2[2] ) return(t2[2]-t1[2]);
266 t1 += t1[1]; t2 += t2[1];
267 }
268 else { /* must be type SYMBOL */
269 t3 = t1 + t1[1]; t4 = t2 + t2[1];
270 t1 += 2; t2 += 2;
271 while ( t1 < t3 && t2 < t4 ) {
272 if ( *t1 != *t2 ) return(*t2-*t1);
273 if ( t1[1] != t2[1] ) return(t2[1]-t1[1]);
274 t1 += 2; t2 += 2;
275 }
276 if ( t1 < t3 ) return(-1);
277 if ( t2 < t4 ) return(1);
278 }
279 }
280 else return(*t2-*t1);
281 }
282 if ( t1 < tstop1 ) return(-1);
283 if ( t2 < tstop2 ) return(1);
284 return(0);
285 }
286
287 #endif
288
289 /*
290 #] ComparePoly :
291 #[ ConvertToPoly :
292 */
293 /**
294 * Converts a generic term to polynomial notation in which there are
295 * only symbols and brackets.
296 * During conversion there will be only symbols. Brackets are stripped.
297 * Objects that need 'translation' are put inside a special compiler
298 * buffer and represented by a symbol. The numbering of the extra
299 * symbols is down from the maximum. In principle there can be a
300 * problem when running into the already assigned ones.
301 * The output overwrites the input.
302 * comlist is the compiler code. Used for the various options
303 */
304
305 static int FirstWarnConvertToPoly = 1;
306
ConvertToPoly(PHEAD WORD * term,WORD * outterm,WORD * comlist,WORD par)307 int ConvertToPoly(PHEAD WORD *term, WORD *outterm, WORD *comlist, WORD par)
308 {
309 WORD *tout, *tstop, ncoef, *t, *r, *tt, *ttwo = 0;
310 int i, action = 0;
311 tt = term + *term;
312 ncoef = ABS(tt[-1]);
313 tstop = tt - ncoef;
314 tout = outterm+1;
315 t = term + 1;
316 if ( comlist[2] == DOALL ) {
317 while ( t < tstop ) {
318 if ( *t == SYMBOL ) {
319 r = t+2;
320 t += t[1];
321 while ( r < t ) {
322 if ( r[1] > 0 ) {
323 *tout++ = SYMBOL;
324 *tout++ = 4;
325 *tout++ = r[0];
326 *tout++ = r[1];
327 }
328 else {
329 tout[1] = SYMBOL;
330 tout[2] = 4;
331 tout[3] = r[0];
332 tout[4] = -1;
333 i = FindSubterm(tout+1);
334 *tout++ = SYMBOL;
335 *tout++ = 4;
336 *tout++ = MAXVARIABLES-i;
337 *tout++ = -r[1];
338 action = 1;
339 }
340 r += 2;
341 }
342 }
343 else if ( *t == DOTPRODUCT ) {
344 r = t + 2;
345 t += t[1];
346 while ( r < t ) {
347 tout[1] = DOTPRODUCT;
348 tout[2] = 5;
349 tout[3] = r[0];
350 tout[4] = r[1];
351 if ( r[2] < 0 ) {
352 tout[5] = -1;
353 }
354 else {
355 tout[5] = 1;
356 }
357 i = FindSubterm(tout+1);
358 *tout++ = SYMBOL;
359 *tout++ = 4;
360 *tout++ = MAXVARIABLES-i;
361 *tout++ = ABS(r[2]);
362 r += 3;
363 action = 1;
364 }
365 }
366 else if ( *t == VECTOR ) {
367 r = t + 2;
368 t += t[1];
369 while ( r < t ) {
370 tout[1] = VECTOR;
371 tout[2] = 4;
372 tout[3] = r[0];
373 tout[4] = r[1];
374 i = FindSubterm(tout+1);
375 *tout++ = SYMBOL;
376 *tout++ = 4;
377 *tout++ = MAXVARIABLES-i;
378 *tout++ = 1;
379 r += 2;
380 action = 1;
381 }
382 }
383 else if ( *t == INDEX ) {
384 r = t + 2;
385 t += t[1];
386 while ( r < t ) {
387 tout[1] = INDEX;
388 tout[2] = 3;
389 tout[3] = r[0];
390 i = FindSubterm(tout+1);
391 *tout++ = SYMBOL;
392 *tout++ = 4;
393 *tout++ = MAXVARIABLES-i;
394 *tout++ = 1;
395 r++;
396 action = 1;
397 }
398 }
399 else if ( *t == HAAKJE) {
400 if ( par ) {
401 tout[0] = 1; tout[1] = 1; tout[2] = 3;
402 *outterm = (tout+3)-outterm;
403 if ( NormPolyTerm(BHEAD outterm) < 0 ) return(-1);
404 tout = outterm + *outterm;
405 tout -= 3;
406 i = t[1]; NCOPY(tout,t,i);
407 ttwo = tout-1;
408 }
409 else { t += t[1]; }
410 }
411 else if ( *t >= FUNCTION ) {
412 i = FindSubterm(t);
413 t += t[1];
414 *tout++ = SYMBOL;
415 *tout++ = 4;
416 *tout++ = MAXVARIABLES-i;
417 *tout++ = 1;
418 action = 1;
419 }
420 else {
421 if ( FirstWarnConvertToPoly ) {
422 MLOCK(ErrorMessageLock);
423 MesPrint("Illegal object in conversion to polynomial notation");
424 MUNLOCK(ErrorMessageLock);
425 FirstWarnConvertToPoly = 0;
426 }
427 return(-1);
428 }
429 }
430 NCOPY(tout,tstop,ncoef)
431 if ( ttwo ) {
432 WORD hh = *ttwo;
433 *ttwo = tout-ttwo;
434 if ( ( i = NormPolyTerm(BHEAD ttwo) ) >= 0 ) i = action;
435 tout = ttwo + *ttwo;
436 *ttwo = hh;
437 *outterm = tout - outterm;
438 }
439 else {
440 *outterm = tout-outterm;
441 if ( ( i = NormPolyTerm(BHEAD outterm) ) >= 0 ) i = action;
442 }
443 }
444 else if ( comlist[2] == ONLYFUNCTIONS ) {
445 while ( t < tstop ) {
446 if ( *t >= FUNCTION ) {
447 if ( comlist[1] == 3 ) {
448 i = FindSubterm(t);
449 t += t[1];
450 *tout++ = SYMBOL;
451 *tout++ = 4;
452 *tout++ = MAXVARIABLES-i;
453 *tout++ = 1;
454 action = 1;
455 }
456 else {
457 for ( i = 3; i < comlist[1]; i++ ) {
458 if ( *t == comlist[i] ) break;
459 }
460 if ( i < comlist[1] ) {
461 i = FindSubterm(t);
462 t += t[1];
463 *tout++ = SYMBOL;
464 *tout++ = 4;
465 *tout++ = MAXVARIABLES-i;
466 *tout++ = 1;
467 action = 1;
468 }
469 else {
470 i = t[1]; NCOPY(tout,t,i);
471 }
472 }
473 }
474 else {
475 i = t[1]; NCOPY(tout,t,i);
476 }
477 }
478 NCOPY(tout,tstop,ncoef)
479 *outterm = tout-outterm;
480 Normalize(BHEAD outterm);
481 i = action;
482 }
483 else {
484 MLOCK(ErrorMessageLock);
485 MesPrint("Illegal internal code in conversion to polynomial notation");
486 MUNLOCK(ErrorMessageLock);
487 i = -1;
488 }
489 return(i);
490 }
491
492 /*
493 #] ConvertToPoly :
494 #[ LocalConvertToPoly :
495 */
496 /**
497 * Converts a generic term to polynomial notation in which there are
498 * only symbols and brackets.
499 * During conversion there will be only symbols. Brackets are stripped.
500 * Objects that need 'translation' are put inside a special compiler
501 * buffer and represented by a symbol. The numbering of the extra
502 * symbols is down from the maximum. In principle there can be a
503 * problem when running into the already assigned ones.
504 * This uses the FindTree for searching in the global tree and
505 * then looks further in the AT.ebufnum. This allows fully parallel
506 * processing. Hence we need no locks. Cannot be used in the same
507 * module as ConvertToPoly.
508 */
509
LocalConvertToPoly(PHEAD WORD * term,WORD * outterm,WORD startebuf,WORD par)510 int LocalConvertToPoly(PHEAD WORD *term, WORD *outterm, WORD startebuf, WORD par)
511 {
512 WORD *tout, *tstop, ncoef, *t, *r, *tt, *ttwo = 0;
513 int i, action = 0;
514 tt = term + *term;
515 ncoef = ABS(tt[-1]);
516 tstop = tt - ncoef;
517 tout = outterm+1;
518 t = term + 1;
519 while ( t < tstop ) {
520 if ( *t == SYMBOL ) {
521 r = t+2;
522 t += t[1];
523 while ( r < t ) {
524 if ( r[1] > 0 ) {
525 *tout++ = SYMBOL;
526 *tout++ = 4;
527 *tout++ = r[0];
528 *tout++ = r[1];
529 }
530 else {
531 tout[1] = SYMBOL;
532 tout[2] = 4;
533 tout[3] = r[0];
534 tout[4] = -1;
535 i = FindLocalSubterm(BHEAD tout+1,startebuf);
536 *tout++ = SYMBOL;
537 *tout++ = 4;
538 *tout++ = MAXVARIABLES-i;
539 *tout++ = -r[1];
540 action = 1;
541 }
542 r += 2;
543 }
544 }
545 else if ( *t == DOTPRODUCT ) {
546 r = t + 2;
547 t += t[1];
548 while ( r < t ) {
549 tout[1] = DOTPRODUCT;
550 tout[2] = 5;
551 tout[3] = r[0];
552 tout[4] = r[1];
553 if ( r[2] < 0 ) {
554 tout[5] = -1;
555 }
556 else {
557 tout[5] = 1;
558 }
559 i = FindLocalSubterm(BHEAD tout+1,startebuf);
560 *tout++ = SYMBOL;
561 *tout++ = 4;
562 *tout++ = MAXVARIABLES-i;
563 *tout++ = ABS(r[2]);
564 r += 3;
565 action = 1;
566 }
567 }
568 else if ( *t == VECTOR ) {
569 r = t + 2;
570 t += t[1];
571 while ( r < t ) {
572 tout[1] = VECTOR;
573 tout[2] = 4;
574 tout[3] = r[0];
575 tout[4] = r[1];
576 i = FindLocalSubterm(BHEAD tout+1,startebuf);
577 *tout++ = SYMBOL;
578 *tout++ = 4;
579 *tout++ = MAXVARIABLES-i;
580 *tout++ = 1;
581 r += 2;
582 action = 1;
583 }
584 }
585 else if ( *t == INDEX ) {
586 r = t + 2;
587 t += t[1];
588 while ( r < t ) {
589 tout[1] = INDEX;
590 tout[2] = 3;
591 tout[3] = r[0];
592 i = FindLocalSubterm(BHEAD tout+1,startebuf);
593 *tout++ = SYMBOL;
594 *tout++ = 4;
595 *tout++ = MAXVARIABLES-i;
596 *tout++ = 1;
597 r++;
598 action = 1;
599 }
600 }
601 else if ( *t == HAAKJE) {
602 if ( par ) {
603 tout[0] = 1; tout[1] = 1; tout[2] = 3;
604 *outterm = (tout+3)-outterm;
605 if ( NormPolyTerm(BHEAD outterm) < 0 ) return(-1);
606 tout = outterm + *outterm;
607 tout -= 3;
608 i = t[1]; NCOPY(tout,t,i);
609 ttwo = tout-1;
610 }
611 else { t += t[1]; }
612 }
613 else if ( *t >= FUNCTION ) {
614 i = FindLocalSubterm(BHEAD t,startebuf);
615 t += t[1];
616 *tout++ = SYMBOL;
617 *tout++ = 4;
618 *tout++ = MAXVARIABLES-i;
619 *tout++ = 1;
620 action = 1;
621 }
622 else {
623 if ( FirstWarnConvertToPoly ) {
624 MLOCK(ErrorMessageLock);
625 MesPrint("Illegal object in conversion to polynomial notation");
626 MUNLOCK(ErrorMessageLock);
627 FirstWarnConvertToPoly = 0;
628 }
629 return(-1);
630 }
631 }
632 NCOPY(tout,tstop,ncoef)
633 if ( ttwo ) {
634 WORD hh = *ttwo;
635 *ttwo = tout-ttwo;
636 if ( ( i = NormPolyTerm(BHEAD ttwo) ) >= 0 ) i = action;
637 tout = ttwo + *ttwo;
638 *ttwo = hh;
639 *outterm = tout - outterm;
640 }
641 else {
642 *outterm = tout-outterm;
643 if ( ( i = NormPolyTerm(BHEAD outterm) ) >= 0 ) i = action;
644 }
645 return(i);
646 }
647
648 /*
649 #] LocalConvertToPoly :
650 #[ ConvertFromPoly :
651
652 Converts a generic term from polynomial notation to the original
653 in which the extra symbols have been replaced by their values.
654 The output is in outterm.
655 We only deal with the extra symbols in the range from < i <= to
656 The output has to be sent to TestSub because it may contain
657 subexpressions when extra symbols have been replaced.
658 */
659
ConvertFromPoly(PHEAD WORD * term,WORD * outterm,WORD from,WORD to,WORD offset,WORD par)660 int ConvertFromPoly(PHEAD WORD *term, WORD *outterm, WORD from, WORD to, WORD offset, WORD par)
661 {
662 WORD *tout, *tstop, *tstop1, ncoef, *t, *r, *tt;
663 int i;
664 /* first = 1; */
665 tt = term + *term;
666 tout = outterm+1;
667 ncoef = ABS(tt[-1]);
668 tstop = tt - ncoef;
669 /*
670 r = t = term + 1;
671 while ( t < tstop ) {
672 if ( *t == SYMBOL ) {
673 tstop1 = t + t[1];
674 tt = t + 2;
675 while ( tt < tstop1 ) {
676 if ( ( *tt < MAXVARIABLES - to )
677 || ( *tt >= MAXVARIABLES - from ) ) {
678 tt += 2;
679 }
680 else break;
681 }
682 if ( tt >= tstop1 ) { t = tstop1; continue; }
683 while ( r < t ) *tout++ = *r++;
684 t += 2;
685 first = 0;
686 while ( t < tstop1 ) {
687 if ( ( *t < MAXVARIABLES - to )
688 || ( *t >= MAXVARIABLES - from ) ) {
689 *tout++ = SYMBOL;
690 *tout++ = 4;
691 *tout++ = *t++;
692 *tout++ = *t++;
693 }
694 else {
695 *tout++ = SUBEXPRESSION;
696 *tout++ = SUBEXPSIZE;
697 *tout++ = MAXVARIABLES - *t++ + offset;
698 *tout++ = *t++;
699 if ( par ) *tout++ = AT.ebufnum;
700 else *tout++ = AM.sbufnum;
701 FILLSUB(tout)
702 }
703 }
704 r = t;
705 }
706 else {
707 t += t[1];
708 }
709 }
710 if ( first ) {
711 i = *term; t = term;
712 NCOPY(outterm,t,i);
713 return(*term);
714 }
715 while ( r < t ) *tout++ = *r++;
716 NCOPY(tout,tstop,ncoef)
717 *outterm = tout-outterm;
718 */
719 t = term + 1;
720 while ( t < tstop ) {
721 if ( *t == SYMBOL ) {
722 tstop1 = t + t[1];
723 tt = t + 2;
724 while ( tt < tstop1 ) {
725 if ( ( *tt < MAXVARIABLES - to )
726 || ( *tt >= MAXVARIABLES - from ) ) {
727 tt += 2;
728 }
729 else {
730 *tout++ = SUBEXPRESSION;
731 *tout++ = SUBEXPSIZE;
732 *tout++ = MAXVARIABLES - *tt++ + offset;
733 *tout++ = *tt++;
734 if ( par ) *tout++ = AT.ebufnum;
735 else *tout++ = AM.sbufnum;
736 FILLSUB(tout)
737 }
738 }
739 r = tout; t += 2;
740 *tout++ = SYMBOL; *tout++ = 0;
741 while ( t < tstop1 ) {
742 if ( ( *t < MAXVARIABLES - to )
743 || ( *t >= MAXVARIABLES - from ) ) {
744 *tout++ = *t++;
745 *tout++ = *t++;
746 }
747 else { t += 2; }
748 }
749 r[1] = tout - r;
750 if ( r[1] <= 2 ) tout = r;
751 }
752 else {
753 i = t[1]; NCOPY(tout,t,i)
754 }
755 }
756 NCOPY(tout,tstop,ncoef)
757 *outterm = tout-outterm;
758 return(*outterm);
759 }
760
761 /*
762 #] ConvertFromPoly :
763 #[ FindSubterm :
764
765 In this routine we look up a variable.
766 If we don't find it we will enter it in the subterm compiler buffer
767 Searching is by tree structure.
768 Adding changes the tree.
769
770 Notice that in TFORM we should be in sequential mode.
771 */
772
FindSubterm(WORD * subterm)773 WORD FindSubterm(WORD *subterm)
774 {
775 WORD old[5], *ss, *term, number;
776 CBUF *C = cbuf + AM.sbufnum;
777 LONG oldCpointer;
778 term = subterm-1;
779 ss = subterm+subterm[1];
780 /*
781 Convert to proper term
782 */
783 old[0] = *term; old[1] = ss[0]; old[2] = ss[1]; old[3] = ss[2]; old[4] = ss[3];
784 ss[0] = 1; ss[1] = 1; ss[2] = 3; ss[3] = 0; *term = subterm[1]+4;
785 /*
786 We may have to add the term to the compiler
787 buffer and then to the tree. This cannot be done in parallel and
788 hence we have to set a lock.
789 */
790 LOCK(AM.sbuflock);
791
792 oldCpointer = C->Pointer-C->Buffer; /* Offset of course !!!!!*/
793 AddRHS(AM.sbufnum,1);
794 AddNtoC(AM.sbufnum,*term,term,8);
795 AddToCB(C,0)
796 /*
797 See whether we have this one already. If not, insert it in the tree.
798 */
799 number = InsTree(AM.sbufnum,C->numrhs);
800 /*
801 Restore old values and return what is needed.
802 */
803 if ( number < (C->numrhs) ) { /* It existed already */
804 C->Pointer = oldCpointer + C->Buffer;
805 C->numrhs--;
806 }
807 else {
808 GETIDENTITY
809 WORD dim = DimensionSubterm(subterm);
810
811 if ( dim == -MAXPOSITIVE ) { /* Give error message but continue */
812 WORD *old = AN.currentTerm;
813 AN.currentTerm = term;
814 MLOCK(ErrorMessageLock);
815 MesPrint("Dimension out of range in %t");
816 MUNLOCK(ErrorMessageLock);
817 AN.currentTerm = old;
818 }
819 /*
820 Store the dimension
821 */
822 C->dimension[number] = dim;
823 }
824 UNLOCK(AM.sbuflock);
825
826 *term = old[0]; ss[0] = old[1]; ss[1] = old[2]; ss[2] = old[3]; ss[3] = old[4];
827 return(number);
828 }
829
830 /*
831 #] FindSubterm :
832 #[ FindLocalSubterm :
833
834 In this routine we look up a variable.
835 If we don't find it we will enter it in the subterm compiler buffer
836 Searching is by tree structure.
837 Adding changes the tree.
838
839 Notice that in TFORM we should be in sequential mode.
840 */
841
FindLocalSubterm(PHEAD WORD * subterm,WORD startebuf)842 WORD FindLocalSubterm(PHEAD WORD *subterm, WORD startebuf)
843 {
844 WORD old[5], *ss, *term, number, i, j, *t1, *t2;
845 CBUF *C = cbuf + AT.ebufnum;
846 term = subterm-1;
847 ss = subterm+subterm[1];
848 /*
849 Convert to proper term
850 */
851 old[0] = *term; old[1] = ss[0]; old[2] = ss[1]; old[3] = ss[2]; old[4] = ss[3];
852 ss[0] = 1; ss[1] = 1; ss[2] = 3; ss[3] = 0; *term = subterm[1]+4;
853 /*
854 First see whether we have this one already in the global buffer.
855 */
856 number = FindTree(AM.sbufnum,term);
857 if ( number > 0 ) goto wearehappy;
858 /*
859 Now look whether it is in the ebufnum between startebuf and numrhs
860 Note however that we need an offset of (numxsymbol-startebuf)
861 */
862 for ( i = startebuf+1; i <= C->numrhs; i++ ) {
863 t1 = C->rhs[i]; t2 = term;
864 if ( *t1 == *t2 ) {
865 j = *t1;
866 while ( *t1 == *t2 && j > 0 ) { t1++; t2++; j--; }
867 if ( j <= 0 ) {
868 number = i-startebuf+numxsymbol;
869 goto wearehappy;
870 }
871 }
872 }
873 /*
874 Now we have to add it to cbuf[AT.ebufnum]
875 */
876 AddRHS(AT.ebufnum,1);
877 AddNtoC(AT.ebufnum,*term,term,9);
878 AddToCB(C,0)
879 number = C->numrhs-startebuf+numxsymbol;
880 wearehappy:
881 *term = old[0]; ss[0] = old[1]; ss[1] = old[2]; ss[2] = old[3]; ss[3] = old[4];
882 return(number);
883 }
884
885 /*
886 #] FindLocalSubterm :
887 #[ PrintSubtermList :
888
889 Prints all the expressions in the subterm compiler buffer.
890 The format is such that they give definitions of the temporary
891 variables of which the contents are stored in this buffer.
892 These variables have the names Z_123 etc.
893 */
894
PrintSubtermList(int from,int to)895 void PrintSubtermList(int from,int to)
896 {
897 UBYTE buffer[80], *out, outbuffer[300];
898 int first, i, ii, inc = 1;
899 WORD *term;
900 CBUF *C = cbuf + AM.sbufnum;
901 /*
902 if ( to < from ) inc = -1;
903 if ( to == from ) inc = 0;
904 */
905 if ( from <= to ) {
906 inc = 1; to += inc;
907 }
908 else {
909 inc = -1; to += inc;
910 }
911 AO.OutFill = AO.OutputLine = outbuffer;
912 AO.OutStop = AO.OutputLine+AC.LineLength;
913 AO.IsBracket = 0;
914 AO.OutSkip = 3;
915
916 if ( AC.OutputMode == FORTRANMODE || AC.OutputMode == PFORTRANMODE ) {
917 TokenToLine((UBYTE *)" ");
918 AO.OutSkip = 7;
919 }
920 else if ( ( AO.Optimize.debugflags & 1 ) == 1 ) {}
921 else if ( AO.OutSkip > 0 ) {
922 for ( i = 0; i < AO.OutSkip; i++ ) TokenToLine((UBYTE *)" ");
923 }
924 i = from;
925 do {
926 if ( ( AO.Optimize.debugflags & 1 ) == 1 ) {
927 TokenToLine((UBYTE *)"id ");
928 for ( ii = 3; ii < AO.OutSkip; ii++ ) TokenToLine((UBYTE *)" ");
929 }
930 /*
931 if ( AC.OutputMode == NORMALFORMAT ) {
932 TokenToLine((UBYTE *)"id ");
933 }
934 */
935 else if ( AC.OutputMode == FORTRANMODE || AC.OutputMode == PFORTRANMODE ) {}
936 else { TokenToLine((UBYTE *)" "); }
937
938 out = StrCopy((UBYTE *)AC.extrasym,buffer);
939 if ( AC.extrasymbols == 0 ) {
940 out = NumCopy(i,out);
941 out = StrCopy((UBYTE *)"_",out);
942 }
943 else if ( AC.extrasymbols == 1 ) {
944 out = AddArrayIndex(i,out);
945 }
946 out = StrCopy((UBYTE *)"=",out);
947 TokenToLine(buffer);
948 term = C->rhs[i];
949 first = 1;
950 if ( *term == 0 ) {
951 out = StrCopy((UBYTE *)"0",buffer);
952 if ( AC.OutputMode != FORTRANMODE && AC.OutputMode != PFORTRANMODE ) {
953 out = StrCopy((UBYTE *)";",out);
954 }
955 TokenToLine(buffer);
956 }
957 else {
958 while ( *term ) {
959 if ( WriteInnerTerm(term,first) ) Terminate(-1);
960 term += *term;
961 first = 0;
962 }
963 if ( AC.OutputMode != FORTRANMODE && AC.OutputMode != PFORTRANMODE ) {
964 out = StrCopy((UBYTE *)";",buffer);
965 TokenToLine(buffer);
966 }
967 }
968 /*
969 There is a problem with FiniLine because it prepares for a
970 continuation line in fortran mode.
971 But the next statement should start on a blank line.
972 */
973 /*
974 FiniLine();
975 if ( AC.OutputMode == FORTRANMODE || AC.OutputMode == PFORTRANMODE ) {
976 AO.OutFill = AO.OutputLine;
977 TokenToLine((UBYTE *)" ");
978 AO.OutSkip = 7;
979 }
980 */
981 if ( AC.OutputMode == FORTRANMODE || AC.OutputMode == PFORTRANMODE ) {
982 AO.OutSkip = 6;
983 FiniLine();
984 AO.OutSkip = 7;
985 }
986 else {
987 FiniLine();
988 }
989 i += inc;
990 } while ( i != to );
991 }
992
993 /*
994 #] PrintSubtermList :
995 #[ PrintExtraSymbol :
996
997 Prints the definition of extra symbol num as the contents
998 of the expression in terms.
999 The parameter par has three options:
1000 EXTRASYMBOL num is interpreted as the number of an extra symbol
1001 REGULARSYMBOL num is interpreted as the number of a symbol.
1002 It could still be an extra symbol.
1003 EXPRESSIONNUMBER num is the number of an expression.
1004 terms contains the rhs expression.
1005 */
1006
PrintExtraSymbol(int num,WORD * terms,int par)1007 void PrintExtraSymbol(int num, WORD *terms,int par)
1008 {
1009 UBYTE buffer[80], *out, outbuffer[300];
1010 int first, i;
1011 WORD *term;
1012
1013 AO.OutFill = AO.OutputLine = outbuffer;
1014 AO.OutStop = AO.OutputLine+AC.LineLength;
1015 AO.IsBracket = 0;
1016
1017 if ( AC.OutputMode == FORTRANMODE || AC.OutputMode == PFORTRANMODE ) {
1018 TokenToLine((UBYTE *)" ");
1019 AO.OutSkip = 7;
1020 }
1021 else if ( ( AO.Optimize.debugflags & 1 ) == 1 ) {
1022 TokenToLine((UBYTE *)"id ");
1023 for ( i = 3; i < AO.OutSkip; i++ ) TokenToLine((UBYTE *)" ");
1024 }
1025 else if ( AO.OutSkip > 0 ) {
1026 for ( i = 0; i < AO.OutSkip; i++ ) TokenToLine((UBYTE *)" ");
1027 }
1028 out = buffer;
1029 switch ( par ) {
1030 case REGULARSYMBOL:
1031 if ( num >= MAXVARIABLES-cbuf[AM.sbufnum].numrhs ) {
1032 num = MAXVARIABLES-num;
1033 }
1034 else {
1035 out = StrCopy(FindSymbol(num),out);
1036 /* out = StrCopy(VARNAME(symbols,num),out); */
1037 break;
1038 }
1039 /* fall through */
1040 case EXTRASYMBOL:
1041 out = StrCopy(FindExtraSymbol(num),out);
1042 /*
1043 out = StrCopy((UBYTE *)AC.extrasym,out);
1044 if ( AC.extrasymbols == 0 ) {
1045 out = NumCopy(num,out);
1046 out = StrCopy((UBYTE *)"_",out);
1047 }
1048 else if ( AC.extrasymbols == 1 ) {
1049 out = AddArrayIndex(num,out);
1050 }
1051 */
1052 break;
1053 case EXPRESSIONNUMBER:
1054 out = StrCopy(EXPRNAME(num),out);
1055 break;
1056 default:
1057 MesPrint("Illegal option in PrintExtraSymbol");
1058 Terminate(-1);
1059 }
1060 out = StrCopy((UBYTE *)"=",out);
1061 TokenToLine(buffer);
1062 term = terms;
1063 first = 1;
1064 if ( *term == 0 ) {
1065 out = StrCopy((UBYTE *)"0",buffer);
1066 TokenToLine(buffer);
1067 }
1068 else {
1069 while ( *term ) {
1070 if ( WriteInnerTerm(term,first) ) Terminate(-1);
1071 term += *term;
1072 first = 0;
1073 }
1074 }
1075 if ( AC.OutputMode != FORTRANMODE && AC.OutputMode != PFORTRANMODE ) {
1076 out = StrCopy((UBYTE *)";",buffer);
1077 TokenToLine(buffer);
1078 }
1079 FiniLine();
1080 }
1081
1082 /*
1083 #] PrintExtraSymbol :
1084 #[ FindSubexpression :
1085
1086 In this routine we look up a subexpression.
1087 If we don't find it we will enter it in the subterm compiler buffer
1088 Searching is by tree structure.
1089 Adding changes the tree.
1090
1091 Notice that in TFORM we should be in sequential mode.
1092 */
1093
FindSubexpression(WORD * subexpr)1094 WORD FindSubexpression(WORD *subexpr)
1095 {
1096 WORD *term, number;
1097 CBUF *C = cbuf + AM.sbufnum;
1098 LONG oldCpointer;
1099
1100 term = subexpr;
1101 while ( *term ) term += *term;
1102 number = term - subexpr;
1103 /*
1104 We may have to add the subexpression to the tree.
1105 This requires a lock.
1106 */
1107 LOCK(AM.sbuflock);
1108
1109 oldCpointer = C->Pointer-C->Buffer; /* Offset of course !!!!!*/
1110 AddRHS(AM.sbufnum,1);
1111 /*
1112 Add the terms to the compiler buffer. Paste on a zero.
1113 */
1114 AddNtoC(AM.sbufnum,number,subexpr,10);
1115 AddToCB(C,0)
1116 /*
1117 See whether we have this one already. If not, insert it in the tree.
1118 */
1119 number = InsTree(AM.sbufnum,C->numrhs);
1120 /*
1121 Restore old values and return what is needed.
1122 */
1123 if ( number < (C->numrhs) ) { /* It existed already */
1124 C->Pointer = oldCpointer + C->Buffer;
1125 C->numrhs--;
1126 }
1127 else {
1128 GETIDENTITY
1129 WORD dim = DimensionExpression(BHEAD subexpr);
1130 /*
1131 Store the dimension
1132 */
1133 C->dimension[number] = dim;
1134 }
1135
1136 UNLOCK(AM.sbuflock);
1137
1138 return(number);
1139 }
1140
1141 /*
1142 #] FindSubexpression :
1143 #[ ExtraSymFun :
1144 */
1145
ExtraSymFun(PHEAD WORD * term,WORD level)1146 int ExtraSymFun(PHEAD WORD *term,WORD level)
1147 {
1148 WORD *oldworkpointer = AT.WorkPointer;
1149 WORD *termout, *t1, *t2, *t3, *tstop, *tend, i;
1150 int retval = 0;
1151 tend = termout = term + *term;
1152 tstop = tend - ABS(tend[-1]);
1153 t3 = t1 = term+1; t2 = termout+1;
1154 /*
1155 First refind the function(s). There is at least one.
1156 */
1157 while ( t1 < tstop ) {
1158 if ( *t1 == EXTRASYMFUN && t1[1] == FUNHEAD+2 ) {
1159 if ( t1[FUNHEAD] == -SNUMBER && t1[FUNHEAD+1] <= numxsymbol
1160 && t1[FUNHEAD+1] > 0 ) {
1161 i = t1[FUNHEAD+1];
1162 }
1163 else if ( t1[FUNHEAD] == -SYMBOL && t1[FUNHEAD+1] < MAXVARIABLES
1164 && t1[FUNHEAD+1] >= MAXVARIABLES-numxsymbol ) {
1165 i = MAXVARIABLES - t1[FUNHEAD+1];
1166 }
1167 else goto nocase;
1168 while ( t3 < t1 ) *t2++ = *t3++;
1169 /*
1170 Now inset the rhs pointer
1171 */
1172 *t2++ = SUBEXPRESSION;
1173 *t2++ = SUBEXPSIZE;
1174 *t2++ = i;
1175 *t2++ = 1;
1176 *t2++ = AM.sbufnum;
1177 FILLSUB(t2)
1178 t3 = t1 = t1 + t1[1];
1179 }
1180 else if ( *t1 == EXTRASYMFUN && t1[1] == FUNHEAD ) {
1181 while ( t3 < t1 ) *t2++ = *t3++;
1182 t3 = t1 = t1 + t1[1];
1183 }
1184 else {
1185 nocase:;
1186 t1 = t1 + t1[1];
1187 }
1188 }
1189 while ( t3 < tend ) *t2++ = *t3++;
1190 *termout = t2 - termout;
1191 AT.WorkPointer = t2;
1192 if ( AT.WorkPointer >= AT.WorkTop ) {
1193 MLOCK(ErrorMessageLock);
1194 MesWork();
1195 MUNLOCK(ErrorMessageLock);
1196 AT.WorkPointer = oldworkpointer;
1197 return(-1);
1198 }
1199 retval = Generator(BHEAD termout,level);
1200 AT.WorkPointer = oldworkpointer;
1201 if ( retval < 0 ) {
1202 MLOCK(ErrorMessageLock);
1203 MesCall("ExtraSymFun");
1204 MUNLOCK(ErrorMessageLock);
1205 }
1206 return(retval);
1207 }
1208
1209 /*
1210 #] ExtraSymFun :
1211 #[ PruneExtraSymbols :
1212 */
1213
PruneExtraSymbols(WORD downto)1214 int PruneExtraSymbols(WORD downto)
1215 {
1216 CBUF *C = cbuf + AM.sbufnum;
1217 if ( downto < C->numrhs && downto >= 0 ) { /* !!!!! */
1218 ClearTree(AM.sbufnum);
1219 C->numrhs = downto;
1220 if ( downto == 0 ) {
1221 C->Pointer = C->Buffer;
1222 }
1223 else {
1224 WORD *w = C->rhs[downto], i;
1225 while ( *w ) w += *w;
1226 C->Pointer = w+1;
1227 for ( i = 1; i <= downto; i++ ) {
1228 InsTree(AM.sbufnum,i);
1229 }
1230 }
1231 }
1232 return(0);
1233 }
1234
1235 /*
1236 #] PruneExtraSymbols :
1237 */
1238