/[ascend]/trunk/ascend/compiler/relation.c
ViewVC logotype

Annotation of /trunk/ascend/compiler/relation.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2689 - (hide annotations) (download) (as text)
Mon Mar 4 11:28:19 2013 UTC (9 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 142455 byte(s)
Fixes bug 567, for the moment. More tests surely required.
1 johnpye 707 /* ASCEND modelling environment
2 jpye 2326 Copyright (C) 2006, 2010 Carnegie Mellon University
3 johnpye 707 Copyright (C) 1990, 1993, 1994 Thomas Guthrie Epperly
4     Copyright (C) 1993, 1994, 1995 Kirk Andre' Abbott
5     Copyright (C) 1996 Benjamin Andrew Allan
6 aw0a 1
7 johnpye 707 This program is free software; you can redistribute it and/or modify
8     it under the terms of the GNU General Public License as published by
9     the Free Software Foundation; either version 2, or (at your option)
10     any later version.
11    
12     This program is distributed in the hope that it will be useful,
13     but WITHOUT ANY WARRANTY; without even the implied warranty of
14     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15     GNU General Public License for more details.
16    
17     You should have received a copy of the GNU General Public License
18 jpye 2648 along with this program. If not, see <http://www.gnu.org/licenses/>.
19 johnpye 707 *//*
20     @file
21     Relation construction routines
22     *//*
23     by Tom Epperly
24     Created: 1/30/90
25     Last in CVS $Revision: 1.32 $ $Date: 1998/03/17 22:09:24 $ $Author: ballan $
26     */
27    
28 aw0a 1 #include <math.h>
29     #include <stdarg.h>
30 jpye 2323 #include <ascend/general/platform.h>
31 jpye 2322 #include <ascend/general/ascMalloc.h>
32 jpye 2323 #include <ascend/general/panic.h>
33 jpye 2018 #include <ascend/general/pool.h>
34     #include <ascend/general/list.h>
35 johnpye 1031 /* #include <general/pairlist.h> */
36 jpye 2018 #include <ascend/general/stack.h>
37     #include <ascend/general/dstring.h>
38 johnpye 1210
39 johnpye 669 #include "expr_types.h"
40 johnpye 399 #include "name.h"
41     #include "nameio.h"
42     #include "instance_enum.h"
43     #include "bintoken.h"
44     #include "exprs.h"
45     #include "exprio.h"
46     #include "value_type.h"
47     #include "evaluate.h"
48     #include "forvars.h"
49     #include "find.h"
50 johnpye 1031 #include "findpath.h"
51 johnpye 399 #include "sets.h"
52     #include "setinstval.h"
53     #include "instance_io.h"
54     #include "extcall.h"
55     #include "relation_util.h"
56     #include "rel_common.h"
57 johnpye 908 #include "rel_blackbox.h"
58 johnpye 399 #include "temp.h"
59     #include "atomvalue.h"
60     #include "mathinst.h"
61     #include "instquery.h"
62     #include "tmpnum.h"
63 johnpye 908 #include "vlist.h"
64 johnpye 399 #include "relation.h"
65 aw0a 1
66     /*
67     * internal form of RelationRelop for lval or rval use.
68     */
69     #define RelRelop(r) ((r)->share->s.relop)
70    
71     #define SUM 1
72     #define PROD 0
73     #ifndef abs
74     #define abs(a) ( ((a)>0) ? (a) : (-(a)) )
75     #endif
76    
77     /*
78     * Some global and exported variables.
79     */
80     struct gl_list_t *g_relation_var_list = NULL;
81    
82     int g_simplify_relations = 1;
83    
84     int g_ExternalNodeStamps=0; /* incremented each time an new external
85     * statement is seen */
86    
87     /* fwd declaration */
88     static union RelationTermUnion
89     *CopyRelationSide(union RelationTermUnion *, unsigned long);
90    
91     #ifdef THIS_IS_AN_UNUSED_FUNCTION
92     static
93     unsigned long ExprLength(register CONST struct Expr *start,
94     register CONST struct Expr *stop)
95     {
96     register unsigned long result=0;
97     while(start!=stop){
98     start = NextExpr(start);
99     result++;
100     }
101     return result;
102     }
103     #endif /* THIS_IS_AN_UNUSED_FUNCTION */
104    
105    
106     static
107     void FigureOutError(struct value_t value,
108     enum relation_errors *err,
109     enum find_errors *ferr)
110     {
111     assert(ValueKind(value)==error_value);
112     *err = find_error;
113     switch(ErrorValue(value)){
114     case type_conflict:
115     case dimension_conflict:
116     case incorrect_name:
117     case incorrect_such_that:
118     case empty_choice:
119     case empty_intersection:
120     case temporary_variable_reused:
121     *ferr = impossible_instance;
122     break;
123     case undefined_value:
124     *ferr = undefined_instance;
125     break;
126     case name_unfound:
127     *ferr = unmade_instance;
128     break;
129     default:
130 johnpye 1064 ASC_PANIC("Unknown error type in FigureOutError.\n");
131 aw0a 1 break;
132     }
133     }
134    
135 johnpye 669 /*-----------------------------------------------------------------------------
136     CREATION AND MANAGEMENT OF RELATION TERMS
137 johnpye 89
138 aw0a 1 It is cheaper to create relation terms in arrays the size of
139     the union than individually because of operating system overhead.
140    
141     Lookout, the tokens have unionized: next they'll want a raise.
142 johnpye 669 */
143    
144 aw0a 1 /*
145     * The define POOL_ALLOCTERM is for people who are pulling terms out
146     * of a pool and promise to return them immediately.
147     */
148    
149     static pool_store_t g_term_pool = NULL;
150     /* A pool_store for 1 expression.
151     * It is expected that objective functions will cause the
152     * largest expressions.
153     * Each time an expression is completed, it will be copied
154     * into an array which can be created already knowing
155     * its proper size. The array will be naturally in postfix.
156     */
157    
158     #define POOL_ALLOCTERM A_TERM(pool_get_element(g_term_pool))
159     /* get a token. Token is the size of the RelationTermUnion. */
160     #ifdef NDEBUG
161     #define PTINIT(x)
162     #else
163     #define PTINIT(x) TermUnionInit(x)
164     #endif
165     #define POOL_RESET pool_clear_store(g_term_pool)
166     /* reset the pool for next expression */
167    
168     #ifndef NDEBUG
169     /*
170     * this function zeros a termunion ptr contents. tu must not be NULL.
171     */
172     static void TermUnionInit(struct relation_term *tu)
173     {
174     memset((char *)tu,0,sizeof(union RelationTermUnion));
175     }
176     #endif
177    
178     static struct {
179     long startcheck;
180     size_t len;
181     size_t cap;
182     struct relation_term **buf;
183     unsigned long *termstack;
184     unsigned long termstackcap;
185     long endcheck;
186     } g_term_ptrs = {1234567890,0,0,NULL,NULL,0,987654321};
187    
188     #define TPBUF_RESET (g_term_ptrs.len=0)
189     /* forget about all the terms in the buffer */
190    
191    
192     /*
193     * Now one can ask why a pool and a buffer both? Couldn't one just
194     * run a big buffer? Well, yes, but how big? Growing a buffer of
195     * complete tokens can cause some system allocators to behave very
196     * poorly. Growing a vector of pointers to tokens is much less
197     * likely to cause the allocator headaches.
198     *
199     * The pool has a good growth mechanism and can handle tokens.
200     * Tradeoff: it is slower to copy the final token data into a
201     * fixed array from pool pointers than from a buffer monolith.
202     */
203     #define TPBUF_INITSIZE 1000
204     /* initial token buffer capacity */
205     #define TPBUF_GROW 1000
206     /* token buffer growth rate */
207    
208     #define RP_LEN 5
209     #if (SIZEOF_VOID_P == 8)
210     #define RP_WID 41
211     #else
212     #define RP_WID 63
213     #endif
214     /* retune rpwid if the size of tokens changes dramatically */
215     #define RP_ELT_SIZE (sizeof(union RelationTermUnion))
216     #define RP_MORE_ELTS 5
217     /* Number of slots filled if more elements needed.
218     So if the pool grows, it grows by RP_MORE_ELTS*RP_WID elements at a time. */
219     #define RP_MORE_BARS 508
220     /* This is the number of pool bar slots to add during expansion.
221     not all the slots will be filled immediately. */
222    
223     /* This function is called at compiler startup time and destroy at shutdown.
224     One could also recall these every time there is a delete all types. */
225     void InitRelInstantiator(void) {
226     if (g_term_pool != NULL || g_term_ptrs.buf != NULL) {
227 johnpye 1064 ASC_PANIC("ERROR: InitRelInstantiator called twice.\n");
228 aw0a 1 }
229     g_term_pool =
230     pool_create_store(RP_LEN, RP_WID, RP_ELT_SIZE, RP_MORE_ELTS, RP_MORE_BARS);
231     if (g_term_pool == NULL) {
232     Asc_Panic(2, "InitRelInstantiator",
233     "ERROR: InitRelInstantiator unable to allocate pool.\n");
234     }
235     g_term_ptrs.buf = (struct relation_term **)
236 johnpye 709 ASC_NEW_ARRAY_CLEAR(union RelationTermUnion *,TPBUF_INITSIZE);
237 aw0a 1 /* don't let the above cast fool you about what's in the array */
238     if (g_term_ptrs.buf == NULL) {
239     Asc_Panic(2, "InitRelInstantiator",
240     "ERROR: InitRelInstantiator unable to allocate memory.\n");
241     }
242     g_term_ptrs.len = 0;
243     g_term_ptrs.cap = TPBUF_INITSIZE;
244     g_term_ptrs.termstackcap = 200;
245 johnpye 709 g_term_ptrs.termstack = ASC_NEW_ARRAY(unsigned long,200);
246 aw0a 1 if (g_term_ptrs.termstack == NULL) {
247     Asc_Panic(2, "InitRelInstantiator",
248     "ERROR: InitRelInstantiator unable to allocate memory.\n");
249     }
250     }
251    
252     /* this function returns NULL when newcap is 0 or when
253     * it is unable to allocate the space requested.
254     */
255     static unsigned long *realloc_term_stack(unsigned long newcap){
256     if (!newcap) {
257     if (g_term_ptrs.termstackcap !=0) {
258     ascfree(g_term_ptrs.termstack);
259     g_term_ptrs.termstack = NULL;
260     g_term_ptrs.termstackcap = 0;
261     }
262     } else {
263     if (newcap >= g_term_ptrs.termstackcap) { /*less than means currently ok */
264     unsigned long *newbuf;
265     newbuf = (unsigned long *)
266     ascrealloc(g_term_ptrs.termstack,(sizeof(unsigned long)*newcap));
267     if (newbuf!=NULL) {
268     g_term_ptrs.termstack = newbuf;
269     g_term_ptrs.termstackcap = newcap;
270     } else {
271     FPRINTF(ASCERR,"Insufficient memory in relation processor\n");
272     return NULL;
273     }
274     }
275     }
276     return g_term_ptrs.termstack;
277     }
278    
279     void DestroyRelInstantiator(void) {
280     assert(g_term_ptrs.buf!=NULL);
281     assert(g_term_pool!=NULL);
282     ascfree(g_term_ptrs.buf);
283     g_term_ptrs.buf = NULL;
284     g_term_ptrs.cap = g_term_ptrs.len = (size_t)0;
285     if (g_term_ptrs.termstackcap != 0) {
286     ascfree(g_term_ptrs.termstack);
287     g_term_ptrs.termstack = NULL;
288     g_term_ptrs.termstackcap = 0;
289     }
290     pool_destroy_store(g_term_pool);
291     g_term_pool = NULL;
292     }
293    
294     void ReportRelInstantiator(FILE *f)
295     {
296     assert(g_term_pool!=NULL);
297     FPRINTF(f,"RelInstantiator ");
298     pool_print_store(f,g_term_pool,0);
299     FPRINTF(f,"RelInstantiator buffer capacity: %lu\n",
300     (unsigned long)g_term_ptrs.cap);
301     }
302    
303     /* The slower expansion process. */
304     static void ExpandTermBuf(struct relation_term *t) {
305     struct relation_term **newbuf;
306     newbuf = (struct relation_term **)ascrealloc(g_term_ptrs.buf,
307     (sizeof(struct relation_term *)*(g_term_ptrs.cap+TPBUF_GROW)));
308     if (newbuf!=NULL) {
309     g_term_ptrs.buf = newbuf;
310     g_term_ptrs.cap += TPBUF_GROW;
311     g_term_ptrs.buf[g_term_ptrs.len] = t;
312     g_term_ptrs.len++;
313     } else {
314     FPRINTF(ASCERR,
315     "ERROR: Relation Instantiator unable to allocate memory.\n");
316     /* we have ignored the term pointer, but somebody else still has it: pool*/
317     }
318     return;
319     }
320    
321     /* Appends term to buffer. if buffer full and can't expand, forgets term.*/
322     static void AppendTermBuf(struct relation_term *t) {
323     if (g_term_ptrs.len < g_term_ptrs.cap) {
324     g_term_ptrs.buf[g_term_ptrs.len++] = t;
325     } else {
326     ExpandTermBuf(t);
327     }
328     return;
329     }
330    
331 johnpye 669 /*------------------------------------------------------------------------------
332     FUNCS TO SIMPLIFY POSTFIX TOKEN LIST
333 aw0a 1
334 johnpye 669 ...before final creation of the token relation array.
335     */
336    
337 aw0a 1 /* returns 1 if term is e_zero, e_real=0.0, or e_int=0 */
338     static int SimplifyTBIsZero(struct relation_term *arg)
339     {
340     if (RelationTermType(arg)==e_real && R_TERM(arg)->value == 0.0) return 1;
341     if (RelationTermType(arg)==e_int && I_TERM(arg)->ivalue == 0) return 1;
342     if (RelationTermType(arg)==e_zero) return 1;
343     return 0;
344     }
345    
346     #ifdef THIS_IS_AN_UNUSED_FUNCTION
347     /* check a termtype, t, for scalarness. return 1 if so, 0 otherwise. */
348     static int SimplifyTBIsScalar(enum Expr_enum t)
349     {
350     return (t <= TOK_SCALAR_HIGH && t >= TOK_SCALAR_LOW);
351     }
352     #endif /* THIS_IS_AN_UNUSED_FUNCTION */
353    
354    
355     /* check a termtype, t, for constantness, return 1 if so, 0 otherwise. */
356     static int SimplifyTBIsConstant(enum Expr_enum t)
357     {
358     return (t <= TOK_CONSTANT_HIGH && t >= TOK_CONSTANT_LOW);
359     }
360    
361     #define ZEROTERM(rtp) SimplifyTBIsZero(rtp)
362     /* check a term pointer, rtp, for scalarness */
363     #define SCALARTERM(t) SimplifyTBIsScalar(t)
364     /* check a termtype, t, for scalarness */
365     #define CONSTANTTERM(t) SimplifyTBIsConstant(t)
366     /* check a termtype, t, for constantness */
367    
368     /*
369     * Attempt to simplify unary functions.
370     * Returns 1 if arg is not constant.
371     * Returns 0 if succeeded, in which case *fn is now morphed to a constant term.
372     * Returns -1 if arg value/dimens are inconsistent with function fn.
373     * Constant arg with numeric value 0 and wild/no dim are coerced quietly
374     * where applicable.
375     *
376     * Cost: O(1).
377     */
378     static int SimplifyTermBuf_Func(struct relation_term *arg,
379     struct relation_term *fn)
380     {
381     CONST dim_type *newdim=NULL;
382     double rval;
383     /* zero constants */
384     if (ZEROTERM(arg)) {
385     switch(FuncId(F_TERM(fn)->fptr)) {
386     case F_LN:
387 johnpye 123 case F_LOG10:
388 aw0a 1 case F_ARCCOSH:
389     /* illegal argument. caller will whine. */
390     return -1;
391     case F_EXP:
392     case F_COSH:
393     if (IsWild(TermDimensions(arg)) ||
394     SameDimen(TermDimensions(arg),Dimensionless())) {
395     arg->t = e_nop;
396     fn->t = e_int;
397     I_TERM(fn)->ivalue = 1;
398     return 0;
399     } else {
400     return -1; /* dimensional incompatibility */
401     }
402     case F_COS:
403     if (IsWild(TermDimensions(arg)) ||
404     SameDimen(TermDimensions(arg),TrigDimension())) {
405     arg->t = e_nop;
406     fn->t = e_int;
407     I_TERM(fn)->ivalue = 1;
408     return 0;
409     } else {
410     return -1; /* dimensional incompatibility */
411     }
412     case F_SIN:
413     case F_TAN:
414     if (IsWild(TermDimensions(arg)) ||
415     SameDimen(TermDimensions(arg),TrigDimension())) {
416     arg->t = e_nop;
417     fn->t = e_int;
418     I_TERM(fn)->ivalue = 0;
419     return 0;
420     } else {
421     return -1; /* dimensional incompatibility */
422     }
423 johnpye 89 #ifdef HAVE_ERF
424 aw0a 1 case F_ERF:
425 johnpye 89 #endif
426 aw0a 1 case F_SINH:
427     case F_ARCSINH:
428     case F_TANH:
429     case F_ARCTANH:
430     if (IsWild(TermDimensions(arg)) ||
431     SameDimen(TermDimensions(arg),Dimensionless())) {
432     arg->t = e_nop;
433     fn->t = e_int;
434     I_TERM(fn)->ivalue = 0; /* dimensionless integer 0 */
435     return 0;
436     } else {
437     return -1; /* dimensional incompatibility */
438     }
439     case F_CUBE:
440     {
441     newdim = CubeDimension(TermDimensions(arg),1);
442     if (newdim != NULL) {
443     arg->t = e_nop;
444     fn->t = e_real;
445     R_TERM(fn)->value = 0.0;
446     R_TERM(fn)->dimensions = newdim;
447     return 0;
448     } else {
449     return -1; /* dimensional incompatibility */
450     }
451     }
452     case F_CBRT:
453     {
454     newdim = ThirdDimension(TermDimensions(arg),1);
455     if (newdim != NULL) {
456     arg->t = e_nop;
457     fn->t = e_real;
458     R_TERM(fn)->value = 0.0;
459     R_TERM(fn)->dimensions = newdim;
460     return 0;
461     } else {
462     return -1; /* dimensional incompatibility */
463     }
464     }
465     case F_SQR:
466     {
467     newdim = SquareDimension(TermDimensions(arg),1);
468     if (newdim != NULL) {
469     arg->t = e_nop;
470     fn->t = e_real;
471     R_TERM(fn)->value = 0.0;
472     R_TERM(fn)->dimensions = newdim;
473     return 0;
474     } else {
475     return -1; /* dimensional incompatibility */
476     }
477     }
478     case F_SQRT:
479     {
480     newdim = HalfDimension(TermDimensions(arg),1);
481     if (newdim != NULL) {
482     arg->t = e_nop;
483     fn->t = e_real;
484     R_TERM(fn)->value = 0.0;
485     R_TERM(fn)->dimensions = newdim;
486     return 0;
487     } else {
488     return -1; /* dimensional incompatibility */
489     }
490     }
491     case F_ARCSIN:
492     case F_ARCTAN:
493     if (IsWild(TermDimensions(arg)) ||
494     SameDimen(TermDimensions(arg),Dimensionless())) {
495     arg->t = e_nop;
496     fn->t = e_real;
497     R_TERM(fn)->value = 0.0;
498     R_TERM(fn)->dimensions = TrigDimension();
499     return 0;
500     } else {
501     return -1; /* dimensional incompatibility */
502     }
503     case F_ARCCOS:
504     if (IsWild(TermDimensions(arg)) ||
505     SameDimen(TermDimensions(arg),Dimensionless())) {
506     arg->t = e_nop;
507     fn->t = e_real;
508     R_TERM(fn)->value = F_PI_HALF;
509     R_TERM(fn)->dimensions = TrigDimension();
510     return 0;
511     } else {
512     return -1; /* dimensional incompatibility */
513     }
514     case F_ABS:
515     case F_HOLD:
516     {
517     newdim = TermDimensions(arg);
518     if (newdim != NULL) {
519     arg->t = e_nop;
520     fn->t = e_real;
521     R_TERM(fn)->value = 0.0;
522     R_TERM(fn)->dimensions = newdim;
523     return 0;
524     } else {
525     return -1; /* dimensional insanity */
526     }
527     }
528     case F_LNM:
529     return 1; /* user could change lnm epsilon. can't simplify. */
530     default:
531     FPRINTF(ASCERR,"Unrecognized function in relation.\n");
532     return 1;
533     }
534     }
535     /* nonzero int or real */
536     if( (arg->t == e_int) || (arg->t == e_real) ) {
537     newdim = NULL;
538     if (arg->t == e_int) {
539     rval = (double)I_TERM(arg)->ivalue;
540     } else {
541     rval = R_TERM(arg)->value;
542     }
543     switch(FuncId(F_TERM(fn)->fptr)) {
544     /* things that take any trig arg, return dimensionless */
545     case F_SIN:
546     case F_COS:
547     case F_TAN:
548     if (IsWild(TermDimensions(arg)) ||
549     SameDimen(TermDimensions(arg),TrigDimension())) {
550     newdim = Dimensionless();
551     } else {
552     return -1; /* dimensional incompatibility */
553     }
554     break; /* go to fixup */
555     /* things that require arg >= 1, return dimless */
556     case F_ARCCOSH:
557     if( rval < 1.0 ) return -1;
558     /* fall through */
559     case F_LN:
560 johnpye 123 case F_LOG10:
561 aw0a 1 if( rval < 0.0 ) return -1;
562     if (IsWild(TermDimensions(arg)) ||
563     SameDimen(TermDimensions(arg),Dimensionless())) {
564     newdim = Dimensionless();
565     } else {
566     return -1; /* dimensional incompatibility */
567     }
568     break; /* go to fixup */
569     /* things that take any exponentiable arg, return dimensionless */
570     case F_EXP:
571     case F_SINH:
572     case F_COSH:
573     if (fabs(rval) > F_LIM_EXP) return -1;
574     /* fall through */
575     /* things that take any arg, return dimensionless */
576     case F_ARCSINH:
577     case F_TANH:
578 johnpye 89 #ifdef HAVE_ERG
579     case F_ERF:
580     #endif
581 aw0a 1 if (IsWild(TermDimensions(arg)) ||
582     SameDimen(TermDimensions(arg),Dimensionless())) {
583     newdim = Dimensionless();
584     } else {
585     return -1; /* dimensional incompatibility */
586     }
587     break;
588     case F_ARCTANH:
589     /* things that take any arg abs <1, return dimensionless */
590     if (fabs(rval) < 1.0 && (IsWild(TermDimensions(arg)) ||
591     SameDimen(TermDimensions(arg),Dimensionless()))) {
592     newdim = Dimensionless();
593     } else {
594     return -1; /* dimensional incompatibility or range */
595     }
596     break;
597     case F_CUBE:
598     {
599     newdim = CubeDimension(TermDimensions(arg),1);
600     if (newdim == NULL || fabs(rval) > F_LIM_CUBE) {
601     return -1; /* dimensional incompatibility */
602     }
603     }
604     break;
605     case F_CBRT:
606     {
607     newdim = ThirdDimension(TermDimensions(arg),1);
608     if (newdim == NULL) {
609     return -1; /* dimensional incompatibility , range*/
610     }
611     break;
612     }
613     case F_SQR:
614     {
615     newdim = SquareDimension(TermDimensions(arg),1);
616     if (newdim == NULL || fabs(rval) > F_LIM_SQR) {
617     return -1; /* dimensional incompatibility , range*/
618     }
619     break;
620     }
621     case F_SQRT:
622     {
623     newdim = HalfDimension(TermDimensions(arg),1);
624     if (newdim == NULL || rval < 0.0) {
625     return -1; /* dimensional incompatibility or range */
626     }
627     break;
628     }
629     /* things that take any trig arg, return dimensionless */
630     case F_ARCSIN:
631     case F_ARCCOS:
632     if ( fabs(rval) <= 1.0 && (IsWild(TermDimensions(arg)) ||
633     SameDimen(TermDimensions(arg),Dimensionless()))) {
634     newdim = TrigDimension();
635     break;
636     } else {
637     return -1; /* dimensional incompatibility */
638     }
639     case F_ARCTAN:
640     if (IsWild(TermDimensions(arg)) ||
641     SameDimen(TermDimensions(arg),Dimensionless())) {
642     newdim = TrigDimension();
643     break;
644     } else {
645     return -1; /* dimensional incompatibility */
646     }
647     case F_ABS:
648     case F_HOLD:
649     newdim = TermDimensions(arg);
650     break;
651     case F_LNM:
652     return 1; /* user could change lnm epsilon. can't simplify. */
653     default:
654     FPRINTF(ASCERR,"Unrecognized function in relation.\n");
655     return 1;
656     }
657     rval = FuncEval(TermFunc(A_TERM(fn)),rval);
658     if (floor(rval)==ceil(rval) && SameDimen(newdim,Dimensionless()) &&
659     fabs(rval) < MAXINTREAL) {
660     fn->t = e_int;
661     I_TERM(fn)->ivalue = (long)floor(rval);
662     } else {
663     fn->t = e_real;
664     R_TERM(fn)->value = rval;
665     R_TERM(fn)->dimensions = newdim;
666     }
667     return 0;
668     }
669     return 1;
670     }
671 johnpye 89
672 aw0a 1 static int ArgsForToken(enum Expr_enum t) {
673     switch (t) {
674     case e_nop:
675     case e_undefined:
676     case e_glassbox:
677     case e_blackbox:
678     case e_opcode:
679     case e_token:
680     case e_zero:
681     case e_real:
682     case e_int:
683     case e_var:
684     return 0;
685     case e_uminus:
686     case e_func:
687     return 1;
688     case e_plus:
689     case e_minus:
690     case e_times:
691     case e_divide:
692     case e_power:
693     case e_ipower:
694     case e_notequal:
695     case e_equal:
696     case e_less:
697     case e_greater:
698     case e_lesseq:
699     case e_greatereq:
700     return 2;
701     case e_maximize:
702     case e_minimize:
703     return 1;
704     default:
705     FPRINTF(ASCERR,"ArgsForToken called with illegal token type.\n");
706     return -1;
707     }
708     }
709    
710 johnpye 669 /**
711 aw0a 1 * first = SimplifyTermBuf_SubExprLimit(ts,b,start,tt)
712     * unsigned long CONST *ts; current term stack
713     * struct relation_term ** CONST b; global term ptr array
714     * unsigned long start; starting index IN STACK ts to find needed args
715     * enum Expr_enum tt; term type of operator you want the subexpr for
716     * long int first; term stack position of rightmost arg outside subexpr
717     *
718     * A little function to find the extent of a postfix subexpression for
719     * the args of an operator term in the termstack/termbuf processing.
720     * Returns -2 if insanity detected. handles nonoperator tt gracefully (-2).
721     *
722     * e.g. cos(v1+v2) * v3
723     * tt = e_times, ts =>
724     * | V1 | V2 | + | cos | V3 | * |
725     * ^--------start = 3
726     * ^--------first = -1
727     *
728     * e.g. v1 * (v2 + v3)
729     * tt = e_plus, ts =>
730     * | V1 | V2 | V3 | + | * |
731     * ^--------start = 2
732     * ^--------first = 0
733     *
734     * O(n) n= subexpr length.
735     */
736     static long
737     SimplifyTermBuf_SubExprLimit(unsigned long CONST *ts,
738     struct relation_term ** CONST buf,
739     unsigned long start,
740     enum Expr_enum tt)
741     {
742     long int first, req_args;
743    
744     first = start;
745     req_args = ArgsForToken(tt);
746     if (first < 0) {
747     FPRINTF(ASCERR,"SimplifyTermBuf_SubExpr given malformed subexpression.\n");
748     }
749    
750     while (first >= 0 && req_args >0) {
751     switch(buf[ts[first]]->t) {
752     case e_zero:
753     case e_real:
754     case e_int:
755     case e_var:
756     req_args--;
757     break;
758     case e_plus:
759     case e_minus:
760     case e_times:
761     case e_divide:
762     case e_power:
763     case e_ipower:
764     req_args++;
765     break;
766     case e_func:
767     case e_uminus:
768     break;
769     default:
770     FPRINTF(ASCERR,
771     "SimplifyTermBuf_SubExpr found illegal argument type (%d).\n",
772     buf[ts[first]]->t);
773     return -2;
774     }
775     first--;
776     }
777     if (first < -1) {
778     FPRINTF(ASCERR,"SimplifyTermBuf_SubExpr found malformed subexpression.\n");
779     }
780     return first;
781     }
782    
783     #ifndef NDEBUG
784     /* some functions to keep assert happy when simplification is in debug */
785     static int check_gt0(unsigned long i) {
786     assert(i);
787     return 1;
788     }
789     static int check_gt1(unsigned long i) {
790     assert(i>1);
791     return 1;
792     }
793     #endif
794 johnpye 89
795 johnpye 669 /**
796 johnpye 89 * A function to simplify the term buffer before copying it into a
797 aw0a 1 * postfix array. Only mandatory dim checking is performed.
798     * Cost: O(n) where n = blen.
799     *
800     * This function is rather large, but simply structured, because speed
801 johnpye 89 * is important.
802 aw0a 1 * This is postfix simplification on the cheap. It could be more aggressive,
803     * but only at potentially quadratic expense.
804     *
805     * int level;
806     * struct relation_term ** CONST b;
807     * CONST unsigned long blen;
808     * They are the original term buffer array and its starting length.
809     * b stays constant, not the data in it!
810     *
811     * (the following level definitions are mostly vapor. see relation.h for true.
812     * level is how far to go in simplification. it is cumulative.
813     * level 0 = do nothing.
814     * level 1 = constant folding
815     * level 2 = zero reductions. A*0 = 0/A =0. A^0=1;
816     * level 3 = converting division by constants into multiplication
817     * level 4 = distributing constants over simple mult. (V*C2)*C1 --> V*C3
818     *
819     * As a side effect, any e_power term that can be resolved to having
820     * an integer exponent is converted to an e_ipower.
821     *
822     * This function is designed to simplifications wrt constants that
823     * are easy to do in postfix. If you want something more clever, you
824     * need to dress up things in infix, simplify, and put back to postfix.
825     * Better you than me, bud.
826     *
827     * At present level > 1 is ignored; we always do 1-3, never 4.
828     *
829     * All this goes on in the termbuf array leaving null pointers behind.
830     * We will compact the array and adjust the length before leaving this
831     * function, so you don't have to care about len changing.
832     * The termbuf pointers are from the pool, so we do not free them
833     * as terms are eliminated.
834     *
835     * Internal doc:
836     * Because C optimizers are pretty damned good, we aren't going to
837     * play pointer games, we will just play subscript of b games.
838     * Note that in flight we create null pointers in the already
839     * visited buffer, but we always have an argument immediately
840     * to the left (b[i-1]) of operator b[i]. If b[i] binary, its
841     * right arg is b[i-1] and its left arg is the first nonnull
842     * entry b[j] to the left of b[i-1] (j<i-1).
843     *
844     * The buffer is in postfix. We have no infix to maintain yet.
845     * Abbreviations in comments:
846     * U - unary operator
847     * B - binary operator
848     * P - any operator
849     * V - e_var arg
850     * A - any arg
851     * C - any constant arg (e_int, e_real)
852     * R - e_real arg
853     * I - e_int arg
854     * N - null pointer
855     * While in flight:
856     | A | A | A | A | A | A | A | termbuf
857     * ^------- top = rightmost we've considered (current).
858     | S | S | S | 0 |
859     * ^----next = next free location to put an index in termstack
860     */
861     static unsigned long SimplifyTermBuf(int level,
862     register struct relation_term ** CONST b,
863     CONST unsigned long blen)
864     {
865     register unsigned long next;
866     register unsigned long *ts; /* term stack, should we need it */
867     unsigned long top;
868     long last;
869     unsigned long right;
870     int early = 0, err;
871     CONST dim_type *newdim;
872     long ival;
873     double rval;
874    
875     if ( level < 1 || !blen ) {
876     realloc_term_stack(0);
877     return blen;
878     }
879     ts = realloc_term_stack(blen);
880     /* stack gets used a lot, so make him locally managed, reusable mem */
881     if (ts==NULL) return blen;
882     /* at any trip through this loop we must be able to guarantee
883     * some simple change, or that the buffer is suitable for
884     * cleanup and return, so that we can handle the rogue operators,
885     * args cleanly.
886     */
887     /* check that stack doesn't start with operator */
888     /* should check that stack doesn't start pos 1 with binary operator */
889     switch (b[0]->t) {
890     case e_var:
891     case e_int:
892     case e_real:
893     case e_zero:
894     break;
895     default:
896     FPRINTF(ASCERR,"Compiler cannot simplify malformed expression\n");
897     return blen;
898     }
899    
900     #ifdef NDEBUG
901 johnpye 89 # define TS_TOP (ts[next-1]) /* term address last pushed */
902     # define TS_LEFT (ts[next-2])
903     /* left hand term address IFF current term is binary and the term at TS_TOP is scalar (not operator) */
904     # define TS_SHIFTPOP ts[next-2] = ts[next-1]; next-- /* overwrite ts_left with ts_top and pop */
905 aw0a 1 #else
906 johnpye 89 # define TS_TOP (check_gt0(next),ts[next-1]) /* term address last pushed */
907     # define TS_LEFT (check_gt1(next),ts[next-2]) /* left hand term address IFF current term is binary and the term at TS_TOP is scalar (not operator) */
908     # define TS_SHIFTPOP assert(next>1); ts[next-2] = ts[next-1]; next-- /* overwrite ts_left with ts_top and pop */
909 aw0a 1 #endif
910     /* keep the above definitions in sync. only difference should be assert. */
911    
912 johnpye 89 #define TS_PUSH(index) ts[next]=(index); next++ /* add a term to the stack */
913     #define TS_POP next-- /* backup the stack */
914     #define TS_POP2 next -= 2 /* backup the stack 2 spots */
915    
916 aw0a 1 for (next=top=0; top < blen; top++) {
917     /* pass through the tokens pointers array */
918     if (b[top]==NULL) continue; /* so we can go through again if we like */
919     /* each case and nested case should be complete in itself for
920     readability. do not use fall throughs */
921     switch (b[top]->t) {
922     case e_var:
923     case e_int:
924     case e_real:
925     case e_zero:
926     TS_PUSH(top);
927     break;
928     case e_nop:
929     b[top] = NULL; /* forget nop */
930     break;
931     case e_func:
932     if ( CONSTANTTERM(b[TS_TOP]->t) ) {
933     /* C U -> C' */
934     if ( (err = SimplifyTermBuf_Func(b[TS_TOP],b[top]) ) != 0 ) {
935     /* not simplified. just push later. whine if needed. */
936     if (err < 0) {
937     FPRINTF(ASCERR,
938     "Can't simplify inconsistent argument to unary function.\n");
939     }
940     } else {
941     b[TS_TOP] = NULL; /* kill old arg, func term was morphed. */
942     TS_POP; /* set up to push morphed func in place of arg */
943     }
944     }
945     TS_PUSH(top); /* for all cases in the ifs */
946     break;
947     case e_uminus:
948     switch (b[TS_TOP]->t) {
949     case e_int:
950     I_TERM(b[TS_TOP])->ivalue = -(I_TERM(b[TS_TOP])->ivalue);
951     b[top] = b[TS_TOP]; /* I - => -I */
952     b[TS_TOP] = NULL;
953     TS_POP;
954     TS_PUSH(top);
955     break;
956     case e_real:
957     R_TERM(b[TS_TOP])->value = -(R_TERM(b[TS_TOP])->value);
958     b[top] = b[TS_TOP]; /* R - => -R */
959     b[TS_TOP] = NULL;
960     TS_POP;
961     TS_PUSH(top);
962     break;
963     case e_zero:
964     b[top] = b[TS_TOP]; /* -0 = 0 */
965     b[TS_TOP] = NULL;
966     TS_POP;
967     TS_PUSH(top);
968     break;
969     default: /* aren't going to distribute or do other funky things */
970     TS_PUSH(top);
971     break;
972     }
973     break;
974 johnpye 89
975 aw0a 1 case e_plus:
976     /* A 0 + => NULL NULL A */
977     if ( ZEROTERM(b[TS_TOP]) ) {
978     /*
979     * Note: we really should be checking the dimens of A to match
980     * with dimens of 0 if e_real, but we are can't yet.
981     */
982     b[top] = b[TS_LEFT]; /* overwrite the + with the A */
983     b[TS_LEFT] = NULL; /* null the A old location */
984     b[TS_TOP] = NULL; /* null old location of 0 */
985     TS_POP2;
986     TS_PUSH(top);
987     break;
988     }
989     switch (b[TS_TOP]->t) {
990     case e_var:
991     if ( ZEROTERM(b[TS_LEFT]) ) {
992     /* 0 V + => NULL NULL V */
993     /*
994     * Note: we really should be checking the dimens of V to match
995     * with dimens of 0 if e_real, but we are don't yet.
996     */
997     b[TS_LEFT] = NULL; /* null the zero term */
998     b[top] = b[TS_TOP]; /* overwrite the + with the V */
999     b[TS_TOP] = NULL; /* null old location of V */
1000     TS_POP2;
1001     TS_PUSH(top);
1002     } else {
1003     TS_PUSH(top);
1004     }
1005     break;
1006     /* 2 constant args? mangle C1 C2 + => C3 of appropriate type,if ok. */
1007     case e_int: /* 0 I +, R I +, I I + */
1008     if ( CONSTANTTERM(b[TS_LEFT]->t) ) {
1009     /* 2 constant args. mangle C2 I1 + => C3 of appropriate type,if ok.*/
1010     if (b[TS_LEFT]->t==e_zero) { /* 0 I + */
1011     b[top] = b[TS_TOP]; /* overwrite the + with the I */
1012     b[TS_LEFT] = NULL; /* null the 0 old location */
1013     b[TS_TOP] = NULL; /* null old location of I */
1014     TS_POP2;
1015     TS_PUSH(top);
1016     break;
1017     }
1018     if (b[TS_LEFT]->t == e_int) { /* I2 I1 + */
1019     I_TERM(b[TS_TOP])->ivalue += I_TERM(b[TS_LEFT])->ivalue;
1020     b[top] = b[TS_TOP]; /* overwrite the + with the I1' */
1021     b[TS_LEFT] = NULL; /* null the I2 old location */
1022     b[TS_TOP] = NULL; /* null old location of I1 */
1023     TS_POP2;
1024     TS_PUSH(top);
1025     break;
1026     }
1027     if ( b[TS_LEFT]->t==e_real &&
1028     ( SameDimen(R_TERM(b[TS_LEFT])->dimensions,Dimensionless()) ||
1029     (IsWild(R_TERM(b[TS_LEFT])->dimensions) &&
1030     R_TERM(b[TS_LEFT])->value == 0.0)
1031     )
1032     ) { /* Ri2(possibly wild 0.0) I1 + => I1' */
1033     if (floor(R_TERM(b[TS_LEFT])->value) ==
1034     ceil(R_TERM(b[TS_LEFT])->value) &&
1035     fabs(R_TERM(b[TS_LEFT])->value) < MAXINTREAL) {
1036     I_TERM(b[TS_TOP])->ivalue +=
1037     (long)floor(R_TERM(b[TS_LEFT])->value);
1038     b[top] = b[TS_TOP]; /* overwrite the + with the I1' */
1039     b[TS_LEFT] = NULL; /* null the R2 old location */
1040     b[TS_TOP] = NULL; /* null old location of I1 */
1041     TS_POP2;
1042     TS_PUSH(top);
1043     break;
1044     } else { /* morph + to R' */
1045     b[top]->t = e_real;
1046     R_TERM(b[top])->dimensions = Dimensionless();
1047     R_TERM(b[top])->value =
1048     R_TERM(b[TS_LEFT])->value + (double)I_TERM(b[TS_TOP])->ivalue;
1049     b[TS_LEFT] = NULL; /* null the R2 old location */
1050     b[TS_TOP] = NULL; /* null old location of I1 */
1051     TS_POP2;
1052     TS_PUSH(top);
1053     break;
1054     }
1055     } else { /* dimensional conflict can't be fixed */
1056     FPRINTF(ASCERR,
1057     "Can't simplify dimensionally inconsistent arguments to +.\n");
1058     TS_PUSH(top);
1059     }
1060     break;
1061     } else { /* non C TS_LEFT */
1062     TS_PUSH(top);
1063     }
1064     break;
1065     case e_real: /* 0 R +, R R +, I R + */
1066     if ( CONSTANTTERM(b[TS_LEFT]->t) ) {
1067     /* 2 constant args. mangle C2 R1 + => C3 of appropriate type,if ok.*/
1068     newdim = CheckDimensionsMatch(TermDimensions(b[TS_TOP]),
1069     TermDimensions(b[TS_LEFT]));
1070     if (newdim == NULL) {
1071     FPRINTF(ASCERR,
1072     "Can't simplify dimensionally inconsistent arguments to +.\n");
1073     TS_PUSH(top);
1074     break;
1075     }
1076     if (b[TS_LEFT]->t==e_zero) { /* 0 R + */
1077     b[top] = b[TS_TOP]; /* overwrite the + with the R */
1078     b[TS_LEFT] = NULL; /* null the 0 old location */
1079     b[TS_TOP] = NULL; /* null old location of R */
1080     TS_POP2;
1081     TS_PUSH(top);
1082     /* if R was wild, it stays wild */
1083     break;
1084     }
1085     if (b[TS_LEFT]->t == e_int) { /* I2 R1 + */
1086     R_TERM(b[TS_TOP])->value += (double)I_TERM(b[TS_LEFT])->ivalue;
1087     R_TERM(b[TS_TOP])->dimensions = newdim;
1088     b[top] = b[TS_TOP]; /* overwrite the + with the R1' */
1089     b[TS_LEFT] = NULL; /* null the I2 old location */
1090     b[TS_TOP] = NULL; /* null old location of R1 */
1091     TS_POP2;
1092     TS_PUSH(top);
1093     /* if R wild, R becomes dimensionless */
1094     break;
1095     }
1096     if ( b[TS_LEFT]->t==e_real ) { /* R2 R1 + => R1', if R1' whole->I1'*/
1097     b[top]->t = e_real;
1098     R_TERM(b[top])->dimensions = newdim;
1099     R_TERM(b[top])->value =
1100     R_TERM(b[TS_LEFT])->value + R_TERM(b[TS_TOP])->value;
1101     b[TS_LEFT] = NULL; /* null the R2 old location */
1102     b[TS_TOP] = NULL; /* null old location of R1 */
1103     TS_POP2;
1104     TS_PUSH(top);
1105     /* if integer valued dimless real, convert to int */
1106     if (floor(R_TERM(b[top])->value) == ceil(R_TERM(b[top])->value)
1107     && SameDimen(R_TERM(b[top])->dimensions,Dimensionless()) &&
1108     fabs(R_TERM(b[top])->value) < MAXINTREAL) {
1109     I_TERM(b[top])->ivalue = (long)R_TERM(b[top])->value;
1110     b[top]->t = e_int;
1111     }
1112     break;
1113     } else { /* dimensional conflict can't be fixed */
1114     FPRINTF(ASCERR,
1115     "Can't simplify dimensionally inconsistent arguments to +.\n");
1116     TS_PUSH(top);
1117     }
1118     break;
1119     } else { /* non C TS_LEFT */
1120     TS_PUSH(top);
1121     }
1122     break; /* end eplus, right arg is e_real */
1123     default: /* tstop is not 0, R, I, V */
1124     TS_PUSH(top);
1125     break;
1126     } /* end argtype switch of e_plus */
1127     break;
1128 johnpye 89
1129 aw0a 1 case e_minus:
1130     /* A 0 - => NULL NULL A */
1131     if ( ZEROTERM(b[TS_TOP]) ) {
1132     /*
1133     * Note: we really should be checking the dimens of A to match
1134     * with dimens of 0 if e_real, but we are can't yet.
1135     */
1136     b[top] = b[TS_LEFT]; /* overwrite the - with the A */
1137     b[TS_LEFT] = NULL; /* null the A old location */
1138     b[TS_TOP] = NULL; /* null old location of 0 */
1139     TS_POP2;
1140     TS_PUSH(top);
1141     break;
1142     }
1143     switch (b[TS_TOP]->t) {
1144     case e_var:
1145     if ( ZEROTERM(b[TS_LEFT]) ) {
1146     /* 0 V - => NULL V uminus */
1147     /*
1148     * Note: we really should be checking the dimens of V to match
1149     * with dimens of 0 if e_real, but we are don't yet.
1150     */
1151     b[TS_LEFT] = NULL; /* null the zero term */
1152     b[top]->t = e_uminus; /* morph - to uminus */
1153     TS_SHIFTPOP; /* reduce 0 V => V */
1154     TS_PUSH(top);
1155     } else {
1156     TS_PUSH(top);
1157     }
1158     break;
1159     /* 2 constant args? mangle C1 C2 - => C3 of appropriate type,if ok. */
1160     case e_int: /* 0 I -, R I -, I I - */
1161     if ( CONSTANTTERM(b[TS_LEFT]->t) ) {
1162     /* 2 constant args. mangle C2 I1 - => C3 of appropriate type,if ok.*/
1163     if (b[TS_LEFT]->t==e_zero) { /* 0 I - */
1164     b[top] = b[TS_TOP]; /* overwrite the - with -I */
1165     I_TERM(b[top])->ivalue = -(I_TERM(b[top])->ivalue);
1166     b[TS_LEFT] = NULL; /* null the 0 old location */
1167     b[TS_TOP] = NULL; /* null old location of I */
1168     TS_POP2;
1169     TS_PUSH(top);
1170     break;
1171     }
1172     if (b[TS_LEFT]->t == e_int) { /* I2 I1 - */
1173     I_TERM(b[TS_TOP])->ivalue =
1174     I_TERM(b[TS_LEFT])->ivalue - I_TERM(b[TS_TOP])->ivalue;
1175     b[top] = b[TS_TOP]; /* overwrite the - with the I1' */
1176     b[TS_LEFT] = NULL; /* null the I2 old location */
1177     b[TS_TOP] = NULL; /* null old location of I1 */
1178     TS_POP2;
1179     TS_PUSH(top);
1180     break;
1181     }
1182     if ( b[TS_LEFT]->t==e_real &&
1183     ( SameDimen(R_TERM(b[TS_LEFT])->dimensions,Dimensionless()) ||
1184     (IsWild(R_TERM(b[TS_LEFT])->dimensions) &&
1185     R_TERM(b[TS_LEFT])->value == 0.0)
1186     )
1187     ) { /* Ri2(possibly wild 0.0) I1 - => I1' */
1188     if (floor(R_TERM(b[TS_LEFT])->value) ==
1189     ceil(R_TERM(b[TS_LEFT])->value) &&
1190     fabs(R_TERM(b[TS_LEFT])->value) < MAXINTREAL) {
1191     I_TERM(b[TS_TOP])->ivalue =
1192     (long)floor(R_TERM(b[TS_LEFT])->value)
1193     - I_TERM(b[TS_TOP])->ivalue;
1194     b[top] = b[TS_TOP]; /* overwrite the + with the I1' */
1195     b[TS_LEFT] = NULL; /* null the R2 old location */
1196     b[TS_TOP] = NULL; /* null old location of I1 */
1197     TS_POP2;
1198     TS_PUSH(top);
1199     break;
1200     } else { /* morph - to R' */
1201     b[top]->t = e_real;
1202     R_TERM(b[top])->dimensions = Dimensionless();
1203     R_TERM(b[top])->value =
1204     R_TERM(b[TS_LEFT])->value - (double)I_TERM(b[TS_TOP])->ivalue;
1205     b[TS_LEFT] = NULL; /* null the R2 old location */
1206     b[TS_TOP] = NULL; /* null old location of I1 */
1207     TS_POP2;
1208     TS_PUSH(top);
1209     break;
1210     }
1211     } else { /* dimensional conflict can't be fixed */
1212     FPRINTF(ASCERR,
1213     "Can't simplify dimensionally inconsistent arguments to -.\n");
1214     TS_PUSH(top);
1215     }
1216     break;
1217     } else { /* non C TS_LEFT */
1218     TS_PUSH(top);
1219     }
1220     break;
1221 johnpye 89
1222 aw0a 1 case e_real: /* 0 R -, R R -, I R - */
1223     if ( CONSTANTTERM(b[TS_LEFT]->t) ) {
1224     /* 2 constant args. mangle C2 R1 - => C3 of appropriate type,if ok.*/
1225     newdim = CheckDimensionsMatch(TermDimensions(b[TS_TOP]),
1226     TermDimensions(b[TS_LEFT]));
1227     if (newdim == NULL) {
1228     FPRINTF(ASCERR,
1229     "Can't simplify dimensionally inconsistent arguments to -.\n");
1230     TS_PUSH(top);
1231     break;
1232     }
1233     if (b[TS_LEFT]->t==e_zero) { /* 0 R - */
1234     b[top] = b[TS_TOP]; /* overwrite the - with the R */
1235     R_TERM(b[top])->value = -(R_TERM(b[top])->value);
1236     b[TS_LEFT] = NULL; /* null the 0 old location */
1237     b[TS_TOP] = NULL; /* null old location of R */
1238     TS_POP2;
1239     TS_PUSH(top);
1240     /* if R was wild, it stays wild */
1241     break;
1242     }
1243     if (b[TS_LEFT]->t == e_int) { /* I2 R1 - */
1244     R_TERM(b[TS_TOP])->value =
1245     (double)I_TERM(b[TS_LEFT])->ivalue - R_TERM(b[TS_TOP])->value;
1246     R_TERM(b[TS_TOP])->dimensions = newdim;
1247     b[top] = b[TS_TOP]; /* overwrite the - with the R1' */
1248     b[TS_LEFT] = NULL; /* null the I2 old location */
1249     b[TS_TOP] = NULL; /* null old location of R1 */
1250     TS_POP2;
1251     TS_PUSH(top);
1252     /* if R wild, R becomes dimensionless */
1253     break;
1254     }
1255     if ( b[TS_LEFT]->t==e_real ) { /* R2 R1 - => R1', if R1' whole->I1'*/
1256     b[top]->t = e_real; /* morph - to R */
1257     R_TERM(b[top])->dimensions = newdim;
1258     R_TERM(b[top])->value =
1259     R_TERM(b[TS_LEFT])->value - R_TERM(b[TS_TOP])->value;
1260     b[TS_LEFT] = NULL; /* null the R2 old location */
1261     b[TS_TOP] = NULL; /* null old location of R1 */
1262     TS_POP2;
1263     TS_PUSH(top);
1264     /* if integer valued dimless real, convert to int */
1265     if (floor(R_TERM(b[top])->value) == ceil(R_TERM(b[top])->value)
1266     && SameDimen(R_TERM(b[top])->dimensions,Dimensionless())
1267     && fabs(R_TERM(b[top])->value) < MAXINTREAL) {
1268     I_TERM(b[top])->ivalue = (long)R_TERM(b[top])->value;
1269     b[top]->t = e_int;
1270     }
1271     break;
1272     } else { /* dimensional conflict can't be fixed */
1273     FPRINTF(ASCERR,
1274     "Can't simplify dimensionally inconsistent arguments to -.\n");
1275     TS_PUSH(top);
1276     }
1277     break;
1278     } else { /* non C TS_LEFT */
1279     TS_PUSH(top);
1280     }
1281     break; /* end eminus, right arg is e_real */
1282     default: /* tstop is not 0, R, I, V */
1283     TS_PUSH(top);
1284     break;
1285     } /* end argtype switch of e_minus */
1286     break;
1287 johnpye 89
1288 aw0a 1 case e_times:
1289     /* needs completing. only C*C done at present. need A*0 reductions */
1290     if ( !CONSTANTTERM(b[TS_LEFT]->t) && !CONSTANTTERM(b[TS_TOP]->t) ) {
1291     /* no constants in sight, go on fast. */
1292     TS_PUSH(top);
1293     break;
1294     } else {
1295     /* some constants in sight, try things. */
1296     if (b[TS_LEFT]->t == e_zero || b[TS_TOP]->t == e_zero) {
1297     /* end of A 0 * and 0 A * => 0 */
1298     ival = SimplifyTermBuf_SubExprLimit(ts,b,next-1,e_times);
1299     if ( ival > -2 ) {
1300     for (last = next-1; last > ival; last--) {
1301     b[ts[last]] = NULL; /* kill the subexpression tokens */
1302     }
1303     next = ival + 1; /* big stack pop */
1304     b[top]->t = e_zero;
1305     R_TERM(b[top])->dimensions = WildDimension();
1306     R_TERM(b[top])->value = 0.0;
1307     TS_PUSH(top);
1308     break;
1309     } else {
1310     /* we had an error in subexpression limit search */
1311     TS_PUSH(top);
1312     break;
1313     }
1314     } /* end of A 0 * and 0 A * */
1315     /* NOTE: here we should be watching for 0.0 e_real and 0 e_int,
1316     * but as yet we don't have the dimen derivation utility to
1317     * check these cases and derive a properly dimensioned e_real 0.
1318     * We are not going to do a dimensionally incorrect shortcut
1319     * implementation. BAA 3/96
1320     */
1321     if ( CONSTANTTERM(b[TS_LEFT]->t) ) { /* C A * =?=> ?*/
1322     /* LEFT is now ereal or e_int because it passed the 0 and C tests */
1323     if ( b[TS_TOP]->t == e_real) { /* C R * => C */
1324     if ( b[TS_LEFT]->t == e_real ) { /* R R * => R */
1325     newdim = SumDimensions(R_TERM(b[TS_TOP])->dimensions,
1326     R_TERM(b[TS_LEFT])->dimensions,1);
1327     if ( newdim == NULL || IsWild(newdim) ) { /* bad dim */
1328     FPRINTF(ASCERR,
1329     "Mult. by wild or fractional dimension constant not folded.\n");
1330     TS_PUSH(top);
1331     break;
1332     } else { /* dim ok. morph etimes to be result. */
1333     rval = R_TERM(b[TS_TOP])->value * R_TERM(b[TS_LEFT])->value;
1334     /* god help us if this overflows... */
1335     b[top]->t = e_real;
1336     R_TERM(b[top])->dimensions = newdim;
1337     R_TERM(b[top])->value = rval;
1338     b[TS_TOP] = NULL;
1339     b[TS_LEFT] = NULL;
1340     TS_POP2;
1341     TS_PUSH(top);
1342     break;
1343     }
1344     } else { /* I R * => R */
1345     rval =
1346     R_TERM(b[TS_TOP])->value * (double)I_TERM(b[TS_LEFT])->ivalue;
1347     /* god help us if this overflows... */
1348     b[top]->t = e_real;
1349     R_TERM(b[top])->dimensions = R_TERM(b[TS_TOP])->dimensions;
1350     R_TERM(b[top])->value = rval;
1351     b[TS_TOP] = NULL;
1352     b[TS_LEFT] = NULL;
1353     TS_POP2;
1354     TS_PUSH(top);
1355     break;
1356     }
1357     #ifndef NDEBUG
1358     FPRINTF(ASCERR,"Unexpected error in Simplification (1).\n");
1359     /* NOT REACHED */
1360     break;
1361     #endif
1362     }
1363     if ( b[TS_TOP]->t == e_int) { /* C I * => C */
1364     if ( b[TS_LEFT]->t == e_real ) { /* R I * => R */
1365     rval =
1366     R_TERM(b[TS_LEFT])->value * (double)I_TERM(b[TS_TOP])->ivalue;
1367     /* god help us if this overflows... */
1368     b[top]->t = e_real;
1369     R_TERM(b[top])->dimensions = R_TERM(b[TS_LEFT])->dimensions;
1370     R_TERM(b[top])->value = rval;
1371     b[TS_TOP] = NULL;
1372     b[TS_LEFT] = NULL;
1373     TS_POP2;
1374     TS_PUSH(top);
1375     break;
1376     } else { /* I I * => I */
1377     rval = (double)I_TERM(b[TS_TOP])->ivalue *
1378     (double)I_TERM(b[TS_LEFT])->ivalue;
1379     if (fabs(rval) < (double)(LONG_MAX/2)) {/*result safely integer*/
1380     b[top]->t = e_int;
1381     I_TERM(b[top])->ivalue =
1382     I_TERM(b[TS_TOP])->ivalue * I_TERM(b[TS_LEFT])->ivalue;
1383     b[TS_TOP] = NULL;
1384     b[TS_LEFT] = NULL;
1385     TS_POP2;
1386     TS_PUSH(top);
1387     break;
1388     } else {
1389     b[top]->t = e_real;
1390     R_TERM(b[top])->dimensions = Dimensionless();
1391     R_TERM(b[top])->value = rval;
1392     b[TS_TOP] = NULL;
1393     b[TS_LEFT] = NULL;
1394     TS_POP2;
1395     TS_PUSH(top);
1396     break;
1397     }
1398     }
1399     #ifndef NDEBUG
1400     FPRINTF(ASCERR,"Unexpected error in Simplification (2).\n");
1401     /* NOT REACHED */
1402     break;
1403     #endif
1404     }
1405     } /* end all C A * alternatives */
1406     /* if here, no simplifications worked,
1407     * though constants exist.
1408     */
1409     TS_PUSH(top);
1410     break;
1411     } /* end of case e_times outermost if */
1412     #ifndef NDEBUG
1413     FPRINTF(ASCERR,"Unexpected error in Simplification (3).\n");
1414     /* NOT REACHED */
1415     break;
1416     #endif
1417 johnpye 89
1418 aw0a 1 case e_divide: /* note: A1 A2 / postfix => A1/A2 infix */
1419     /* needs completing only does C/C at present. needs to do 0/A. */
1420     if ( !CONSTANTTERM(b[TS_LEFT]->t) && !CONSTANTTERM(b[TS_TOP]->t) ) {
1421     /* no constants in sight, go on fast. */
1422     TS_PUSH(top);
1423     break;
1424     } else {
1425     /* some constants in sight, try things. */
1426     if (b[TS_LEFT]->t == e_zero && b[TS_TOP]->t != e_zero) {
1427     /* 0 A / => 0 */
1428     ival = SimplifyTermBuf_SubExprLimit(ts,b,next-1,e_divide);
1429     if ( ival > -2 ) {
1430     for (last = next-1; last > ival; last--) {
1431     b[ts[last]] = NULL; /* kill the subexpression tokens */
1432     }
1433     next = ival + 1; /* big stack pop, could be pop2 */
1434     b[top]->t = e_zero;
1435     R_TERM(b[top])->dimensions = WildDimension();
1436     R_TERM(b[top])->value = 0.0;
1437     TS_PUSH(top);
1438     break;
1439     } else {
1440     /* we had an error in subexpression limit search */
1441     TS_PUSH(top);
1442     break;
1443     }
1444     } /* end of 0 A / */
1445     /* NOTE: here we should be watching for 0.0 e_real and 0 e_int,
1446     * but as yet we don't
1447     * check these cases and derive a properly dimensioned e_real 0.
1448     * We are not going to do a dimensionally incorrect shortcut
1449     * implementation. BAA 3/96
1450     */
1451     if ( ZEROTERM(b[TS_TOP]) ) {
1452     /* trap A/0 out */
1453     FPRINTF(ASCERR,"Division by constant 0 not simplified.\n");
1454     top = blen;
1455     early = 1; /* set flag that we punted. */
1456     TS_PUSH(top);
1457     break;
1458     } /* end of A/0 out */
1459     if ( CONSTANTTERM(b[TS_LEFT]->t) ) { /* C A / =?=> ?*/
1460     /* LEFT is now R or I because it passed the 0 and C tests */
1461     if ( b[TS_TOP]->t == e_real) { /* C R / => C */
1462     if ( b[TS_LEFT]->t == e_real ) { /* R R / => R */
1463     newdim = DiffDimensions(R_TERM(b[TS_LEFT])->dimensions,
1464     R_TERM(b[TS_TOP])->dimensions,1);
1465     if ( newdim == NULL || IsWild(newdim) ) { /* bad dim */
1466     FPRINTF(ASCERR,
1467     "Div. by wild or fractional dimension constant not folded.\n");
1468     TS_PUSH(top);
1469     break;
1470     } else { /* dim ok. morph edivide to be result. */
1471     rval = R_TERM(b[TS_LEFT])->value / R_TERM(b[TS_TOP])->value;
1472     /* god help us if this overflows/underflows... */
1473     b[top]->t = e_real;
1474     R_TERM(b[top])->dimensions = newdim;
1475     R_TERM(b[top])->value = rval;
1476     b[TS_TOP] = NULL;
1477     b[TS_LEFT] = NULL;
1478     TS_POP2;
1479     TS_PUSH(top);
1480     break;
1481     }
1482     } else { /* I R / => R */
1483     rval =
1484     ((double)I_TERM(b[TS_LEFT])->ivalue) /R_TERM(b[TS_TOP])->value;
1485     /* god help us if this overflows... */
1486     b[top]->t = e_real;
1487     R_TERM(b[top])->dimensions =
1488     DiffDimensions(Dimensionless(),
1489     R_TERM(b[TS_TOP])->dimensions,0);
1490     /* diff dimens not checked here because top is dimensionless */
1491     R_TERM(b[top])->value = rval;
1492     b[TS_TOP] = NULL;
1493     b[TS_LEFT] = NULL;
1494     TS_POP2;
1495     TS_PUSH(top);
1496     break;
1497     }
1498     #ifndef NDEBUG
1499     FPRINTF(ASCERR,"Unexpected error in Simplification (4).\n");
1500     /* NOT REACHED */
1501     break;
1502     #endif
1503     }
1504     if ( b[TS_TOP]->t == e_int) { /* C I / => C */
1505     if ( b[TS_LEFT]->t == e_real ) { /* R I / => R */
1506     rval =
1507     R_TERM(b[TS_LEFT])->value / (double)I_TERM(b[TS_TOP])->ivalue;
1508     /* god help us if this overflows... */
1509     b[top]->t = e_real;
1510     R_TERM(b[top])->dimensions = R_TERM(b[TS_LEFT])->dimensions;
1511     R_TERM(b[top])->value = rval;
1512     b[TS_TOP] = NULL;
1513     b[TS_LEFT] = NULL;
1514     TS_POP2;
1515     TS_PUSH(top);
1516     break;
1517     } else { /* I I / => R! Integer division is NOT allowed */
1518     rval = (double)I_TERM(b[TS_LEFT])->ivalue;
1519     rval /= (double)I_TERM(b[TS_TOP])->ivalue;
1520     b[top]->t = e_real;
1521     R_TERM(b[top])->dimensions = Dimensionless();
1522     R_TERM(b[top])->value = rval;
1523     b[TS_TOP] = NULL;
1524     b[TS_LEFT] = NULL;
1525     TS_POP2;
1526     TS_PUSH(top);
1527     break;
1528     }
1529     #ifndef NDEBUG
1530     FPRINTF(ASCERR,"Unexpected error in Simplification (5).\n");
1531     /* NOT REACHED */
1532     break;
1533     #endif
1534     }
1535     } /* end all C A / alternatives */
1536     if ( CONSTANTTERM(b[TS_TOP]->t) ) { /* A C / => A (1/C) * */
1537     /* we screened out 0 above, so its int or real */
1538     if (b[TS_TOP]->t == e_real) { /* A R / => A R * */
1539     rval = 1/R_TERM(b[TS_TOP])->value;
1540     /* god help us if this overflows... */
1541     b[top]->t = e_times; /* morph / to * */
1542     /* flip dimens */
1543     R_TERM(b[TS_TOP])->dimensions =
1544     DiffDimensions(Dimensionless(),R_TERM(b[TS_TOP])->dimensions,0);
1545     /* diff dimens not checked here because top is dimensionless */
1546     R_TERM(b[TS_TOP])->value = rval; /* flip value */
1547     TS_PUSH(top);
1548     break;
1549     }
1550     if (b[TS_TOP]->t == e_int) { /* A I / => A I * */
1551     rval = 1.0/(double)I_TERM(b[TS_TOP])->ivalue;
1552     /* god help us if this overflows... */
1553     b[top]->t = e_times; /* morph / to * */
1554     /* flip dimens */
1555     b[TS_TOP]->t = e_real; /* morph int to real */
1556     R_TERM(b[TS_TOP])->value = rval; /* flip value */
1557     R_TERM(b[TS_TOP])->dimensions = Dimensionless();
1558     TS_PUSH(top);
1559     break;
1560     }
1561     } /* end of morphing A/C => A*c' */
1562     /* if here, no simplifications worked,
1563     * though constants exist.
1564     */
1565     TS_PUSH(top);
1566     break;
1567     } /* end of case e_divide outermost if */
1568     /* NOT REACHED */
1569     #ifndef NDEBUG
1570     FPRINTF(ASCERR,"Unexpected error in Simplification (6).\n");
1571     break;
1572     #endif
1573     case e_power: /* DANGER! WILL ROBINSON, DANGER! possible fall through */
1574     /* exponents must be dimensionless to make any sense */
1575     if (b[TS_TOP]->t == e_zero || b[TS_TOP]->t == e_int ||
1576     (b[TS_TOP]->t == e_real &&
1577     ( SameDimen(R_TERM(b[TS_TOP])->dimensions,Dimensionless()) ||
1578     IsWild(R_TERM(b[TS_TOP])->dimensions) ) &&
1579     floor(R_TERM(b[TS_TOP])->value)==ceil(R_TERM(b[TS_TOP])->value) &&
1580     fabs(R_TERM(b[TS_TOP])->value) < MAXINTREAL)
1581     ) { /* big if ipowerable */
1582     if (b[TS_TOP]->t == e_real) { /* morph real to int */
1583     b[TS_TOP]->t = e_int;
1584     I_TERM(b[TS_TOP])->ivalue = (long)R_TERM(b[TS_TOP])->value;
1585     }
1586     /* e_zero and e_int are appropriate to ipower and need no morph */
1587     b[top]->t = e_ipower; /* morph to ipower and fall through */
1588     /* FALL THROUGH! FALL THROUGH! FALL THROUGH! FALL THROUGH! */
1589     /* we aren't supposed to allow fall, but this is really the
1590     most perfect place to do power=>ipower conversion.
1591     Note that very large exponent values may be impossible later. */
1592     } else {
1593     /* still need to code C^R case. A^0 promoted to ipow, not here */
1594     if ( CONSTANTTERM(b[TS_LEFT]->t) && CONSTANTTERM(b[TS_TOP]->t) ) {
1595     /* C is either 0, int, or real. R is nonintegral (or damn big) real.
1596     Because R is not integer, C must be positive and dimensionless,
1597     and also small enough not to overflow: C > 1 =>
1598     check for pow(DBL_MAX,1/R) > R */
1599     if ( !SameDimen(R_TERM(b[TS_TOP])->dimensions,Dimensionless()) &&
1600     !IsWild(R_TERM(b[TS_TOP])->dimensions) ) {
1601     FPRINTF(ASCERR,"Illegal dimensioned exponent found in power.\n");
1602     top=blen;
1603     early = 1; /* set flag that we punted. */
1604     break;
1605     }
1606     if (b[TS_LEFT]->t == e_zero) { /* 0^R, R!=0 => 1 */
1607     b[top]->t = e_int;
1608     I_TERM(b[top])->ivalue = 1;
1609     b[TS_TOP] = NULL;
1610     b[TS_LEFT] = NULL;
1611     TS_POP2;
1612     TS_PUSH(top);
1613     break;
1614     /* end of 0^R */
1615     } else {
1616     if (b[TS_LEFT]->t == e_real) { /* R^R */
1617     if ( !SameDimen(R_TERM(b[TS_LEFT])->dimensions,Dimensionless())
1618     && !IsWild(R_TERM(b[TS_LEFT])->dimensions) ) {
1619     /* can happen on very large exponents */
1620     FPRINTF(ASCERR,
1621     "Illegal dimensioned base raised to nonintegral power.\n");
1622     top = blen;
1623     early = 1; /* set flag that we punted. */
1624     break;
1625     } else { /* R(dimless)^R , err if R ln(C) > ln(DBL_MAX) */
1626     if (R_TERM(b[TS_LEFT])->value < 0) {
1627     /* can happen on very large exponents */
1628     FPRINTF(ASCERR,
1629     "Illegal negative base raised to nonintegral power.\n");
1630     top = blen;
1631     early = 1; /* set flag that we punted. */
1632     break;
1633     }
1634     if (R_TERM(b[TS_LEFT])->value == 0.0) {
1635     /* R!=0, 0^R = 1 */
1636     b[top]->t = e_int;
1637     I_TERM(b[top])->ivalue = 0;
1638     b[TS_TOP] = NULL;
1639     b[TS_LEFT] = NULL;
1640     TS_POP2;
1641     TS_PUSH(top);
1642     break;
1643     }
1644     if ( R_TERM(b[TS_TOP])->value*log(R_TERM(b[TS_LEFT])->value) >
1645     F_LIM_EXP) { /* overflow in result. let solver do it */
1646     TS_PUSH(top);
1647     break;
1648     } else {
1649     b[top]->t = e_real;
1650     R_TERM(b[top])->dimensions = Dimensionless();
1651     R_TERM(b[top])->value =
1652     pow(R_TERM(b[TS_LEFT])->value,R_TERM(b[TS_TOP])->value);
1653     b[TS_TOP] = NULL;
1654     b[TS_LEFT] = NULL;
1655     TS_POP2;
1656     TS_PUSH(top);
1657     break;
1658     }
1659     }
1660     /* end of R^R */
1661     } else { /* I^R */
1662     if (I_TERM(b[TS_LEFT])->ivalue < 0) {
1663     /* can happen on very large exponents */
1664     FPRINTF(ASCERR,
1665     "Illegal negative base raised to nonintegral power.\n");
1666     top = blen;
1667     early = 1; /* set flag that we punted. */
1668     break;
1669     }
1670     if (I_TERM(b[TS_LEFT])->ivalue == 0) {
1671     /* R!=0, 0^R = 1 */
1672     b[top]->t = e_int;
1673     I_TERM(b[top])->ivalue = 0;
1674     b[TS_TOP] = NULL;
1675     b[TS_LEFT] = NULL;
1676     TS_POP2;
1677     TS_PUSH(top);
1678     break;
1679     }
1680     if ( R_TERM(b[TS_TOP])->value *
1681     log((double)I_TERM(b[TS_LEFT])->ivalue) > F_LIM_EXP) {
1682     /* overflow in result. let solver do it */
1683     TS_PUSH(top);
1684     break;
1685     } else {
1686     b[top]->t = e_real;
1687     R_TERM(b[top])->dimensions = Dimensionless();
1688     R_TERM(b[top])->value =
1689     pow((double)I_TERM(b[TS_LEFT])->ivalue,
1690     R_TERM(b[TS_TOP])->value);
1691     b[TS_TOP] = NULL;
1692     b[TS_LEFT] = NULL;
1693     TS_POP2;
1694     TS_PUSH(top);
1695     break;
1696     }
1697     /* end of I^R */
1698     }
1699     /* end of I,R ^R */
1700     }
1701     /* end of 0,I,R ^R */
1702     } else {
1703     TS_PUSH(top);
1704     /* remaining A^A2 where A2 => R or V or expr */
1705     }
1706     /* end of not morphing to ipower */
1707     break;
1708     }
1709     /* FALL THROUGH if morphing to ipower test succeeded */
1710 johnpye 89
1711 aw0a 1 case e_ipower:
1712     if ( ZEROTERM(b[TS_TOP]) ) {
1713     /* A^0 */
1714     if ( ZEROTERM(b[TS_LEFT]) ) {
1715     /* 0^0 */
1716     FPRINTF(ASCERR,"0 raised to 0 power is undefined.\n");
1717     top=blen;
1718     early = 1; /* set flag that we punted. */
1719     break;
1720     } else {
1721     /* A^0 => 1 */
1722     ival = SimplifyTermBuf_SubExprLimit(ts,b,next-1,e_ipower);
1723     if ( ival > -2 ) {
1724     for (last = next-1; last > ival; last--) {
1725     b[ts[last]] = NULL; /* kill the subexpression tokens */
1726     }
1727     next = ival + 1; /* big stack pop */
1728     b[top]->t = e_int;
1729     I_TERM(b[top])->ivalue = 1;
1730     TS_PUSH(top);
1731     break;
1732     } else {
1733     /* we had an error */
1734     TS_PUSH(top);
1735     break;
1736     }
1737     }
1738     } else { /* A^I, I!=0, A!=0 => R or I as needed */
1739     if (CONSTANTTERM(b[TS_LEFT]->t)) { /* C^I */
1740     if (b[TS_LEFT]->t == e_real) { /* R^I */
1741     /* err if I*ln(R) > ln(DBL_MAX) */
1742     if ( I_TERM(b[TS_TOP])->ivalue*log(fabs(R_TERM(b[TS_LEFT])->value))
1743     > F_LIM_EXP) { /* overflow in result. let solver do it */
1744     TS_PUSH(top);
1745     break;
1746     } else {
1747     ival = I_TERM(b[TS_TOP])->ivalue;
1748     newdim = PowDimension(ival,R_TERM(b[TS_LEFT])->dimensions,1);
1749     if (newdim==NULL) {
1750     FPRINTF(ASCERR,
1751     "Dimensional overflow in exponentiation of constant.\n");
1752     TS_PUSH(top);
1753     break;
1754     }
1755     b[top]->t = e_real;
1756     R_TERM(b[top])->dimensions = newdim;
1757     R_TERM(b[top])->value =
1758     asc_ipow(R_TERM(b[TS_LEFT])->value,(int)ival);
1759     /* cast of ival is accurate if newdim was ok */
1760     b[TS_TOP] = NULL;
1761     b[TS_LEFT] = NULL;
1762     TS_POP2;
1763     TS_PUSH(top);
1764     break;
1765     }
1766     /* end of R^I */
1767     } else { /* I^I */
1768     ival = I_TERM(b[TS_TOP])->ivalue;
1769     if ( ival * log((double)abs(I_TERM(b[TS_LEFT])->ivalue))
1770     > F_LIM_EXP) {
1771     /* overflow in result. let solver do it */
1772     TS_PUSH(top);
1773     break;
1774     }
1775     if (abs(ival) < INT_MAX) { /* this could be a little better */
1776     rval = asc_ipow((double)I_TERM(b[TS_LEFT])->ivalue,
1777 ballan 2519 (int)I_TERM(b[TS_TOP])->ivalue);
1778 aw0a 1 if (fabs(rval) > MAXINTREAL || floor(rval)!=ceil(rval) ) {
1779     b[top]->t = e_real;
1780     R_TERM(b[top])->dimensions = Dimensionless();
1781     R_TERM(b[top])->value = rval;
1782     } else { /* can be an int safely */
1783     b[top]->t = e_int;
1784     I_TERM(b[top])->ivalue = (long)rval;
1785     }
1786     b[TS_TOP] = NULL;
1787     b[TS_LEFT] = NULL;
1788     TS_POP2;
1789     TS_PUSH(top);
1790     break;
1791     } else {
1792     /* exponent to damn big */
1793     TS_PUSH(top);
1794     break;
1795     }
1796     /* end of I^I */
1797     } /* end of C^I */
1798     } else {
1799     TS_PUSH(top);
1800     break;
1801     }
1802     #ifndef NDEBUG
1803     FPRINTF(ASCERR,"Unexpected error in Simplification (7).\n");
1804     break; /* NOT REACHED */
1805     #endif
1806     }
1807     #ifndef NDEBUG
1808     FPRINTF(ASCERR,"Unexpected error in Simplification (8).\n");
1809     break; /* NOT REACHED */
1810     #endif
1811     /* end e_ipower */
1812 johnpye 89
1813 aw0a 1 /* all the following are bogus in instantiated tokens at this time. (2/96)
1814     * e_subexpr,e_const,e_par,
1815     * e_card,e_choice,e_sum,e_prod,e_union,e_inter,e_in,e_st,
1816     * e_glassbox,e_blackbox,e_opcode,e_token,
1817     * e_or,e_and,e_boolean,e_set,e_symbol,
1818     * e_equal,e_notequal,e_less,e_greater,e_lesseq,e_greatereq,e_not,
1819     * e_qstring,
1820     * e_maximize,e_minimize,
1821     * e_undefined
1822     */
1823     default:
1824     FPRINTF(ASCERR,"Unexpected token in relation simplification.\n");
1825     FPRINTF(ASCERR,"Returning early.\n");
1826     top=blen;
1827     early = 1; /* flag that we punted. */
1828     break;
1829     }
1830     }
1831     /*
1832     * cleanup reduced expression, if needed.
1833     * after the for loop, next is now the length of the postfix expression,
1834     * or garbage if early is true.
1835     */
1836     if (blen <= next) return blen; /* no simplification, oh well. */
1837     right = 0;
1838     while (right < blen && b[right] != NULL) right++; /* find first null */
1839     for(last = right; right < blen; right++) { /* shift nonnulls left */
1840     if (b[right] != NULL) {
1841     b[last] = b[right];
1842     last++;
1843     }
1844     }
1845 johnpye 89 if (!early && last != (long)next) {
1846 aw0a 1 FPRINTF(ASCERR,"Confusing token counts in Simplify\n");
1847     }
1848     right = last;
1849 johnpye 89 while (last<(long)blen) { /* null remainder, if any, of pointers */
1850 aw0a 1 b[last] = NULL;
1851     last++;
1852     }
1853     return right;
1854     }
1855     /* END SimplifyTermBuf */
1856 johnpye 669
1857    
1858 aw0a 1 struct relation_side_temp {
1859     unsigned long length;
1860     union RelationTermUnion *side;
1861     };
1862    
1863     static struct relation_term
1864     *InfixArr_MakeSide(CONST struct relation_side_temp *, int *);
1865     /* forward declaration */
1866    
1867 johnpye 669 /** returns 1 if converting buf is successful
1868 aw0a 1 * returns 0 if buf empty or insufficient memory.
1869     * The structure tmp given is filled with an array of terms
1870     * and its length. You must free the array if you decide you
1871     * don't want it. We don't care how the structure is initialized.
1872     */
1873     static int ConvertTermBuf(struct relation_side_temp *tmp)
1874     {
1875     union RelationTermUnion *arr = NULL;
1876     unsigned long len,c;
1877    
1878     len = SimplifyTermBuf(g_simplify_relations,g_term_ptrs.buf,g_term_ptrs.len);
1879     if (len < 1) return 0;
1880 johnpye 709 arr = ASC_NEW_ARRAY(union RelationTermUnion,len);
1881 aw0a 1 if (arr==NULL) {
1882     FPRINTF(ASCERR,"Create Token Relation: Insufficient memory :-(.\n");
1883     return 0;
1884     }
1885     for (c=0; c<len; c++) {
1886     arr[c] = *(UNION_TERM(g_term_ptrs.buf[c]));
1887     }
1888     tmp->side = arr;
1889     tmp->length = len;
1890     return 1;
1891     }
1892    
1893 johnpye 669 /**
1894 aw0a 1 * usually we want to reset both simultaneously. reset our
1895     * pooling and buffering data.
1896     */
1897     static
1898     void DestroyTermList(void) {
1899     POOL_RESET;
1900     TPBUF_RESET;
1901     }
1902    
1903 johnpye 669 /**
1904 johnpye 908 create a term from the pool
1905 johnpye 669 */
1906 aw0a 1 static struct relation_term *CreateOpTerm(enum Expr_enum t)
1907     {
1908     struct relation_term *term;
1909     term = POOL_ALLOCTERM;
1910     assert(term!=NULL);
1911     PTINIT(term);
1912     term->t = t;
1913     if (t==e_uminus) {
1914     U_TERM(term)->left = NULL;
1915     } else {
1916     B_TERM(term)->left = NULL;
1917     B_TERM(term)->right = NULL;
1918     }
1919     return term;
1920     }
1921    
1922 johnpye 669 /** create a term from the pool, inserting it
1923 aw0a 1 * in pointer sorted order on g_relation_var_list.
1924     * Note that this and ModifyTokenRelationPointers are the
1925     * only places where the sort
1926     * order of the var list matters.
1927     * In fact, in most cases we could equally afford
1928     * linear search and that would give us repeatability
1929     * across platforms and runs since the vars will be
1930     * then encountered in a constant order determined
1931     * by how the user wrote the equation.
1932     * Needs consideration, especially in light of
1933     * potential to improve relation sharing.
1934     * In particular, we could then easily share
1935     * in a fine-grained manner those relations with
1936     * only a single index involved and no internal sums/products,
1937     * such as f[i] = x[i]*Ftot; in[i].f = out[i].f;
1938     * x = hold(x);
1939     * which could be pretty darn common forms.
1940     */
1941     static struct relation_term *CreateVarTerm(CONST struct Instance *i)
1942     {
1943     struct relation_term *term;
1944     unsigned long pos;
1945 jds 114 if (0 != (pos = gl_search(g_relation_var_list,i,(CmpFunc)CmpP))) {
1946 aw0a 1 /* find var if already on relations var list */
1947     term = POOL_ALLOCTERM;
1948     assert(term!=NULL);
1949     PTINIT(term);
1950     term->t = e_var;
1951     V_TERM(term) -> varnum = pos;
1952     } else {
1953     /* or add it to the var list */
1954     gl_append_ptr(g_relation_var_list,(VOIDPTR)i);
1955     term = POOL_ALLOCTERM;
1956     assert(term!=NULL);
1957     PTINIT(term);
1958     term->t = e_var;
1959     V_TERM(term) -> varnum = gl_length(g_relation_var_list);
1960     }
1961     return term;
1962     }
1963    
1964 johnpye 669 /** create a term from the pool */
1965 aw0a 1 static struct relation_term *CreateIntegerTerm(long int v)
1966     {
1967     struct relation_term *term;
1968     term = POOL_ALLOCTERM;
1969     assert(term!=NULL);
1970     PTINIT(term);
1971     term->t = e_int;
1972     I_TERM(term) -> ivalue = v;
1973     return term;
1974     }
1975    
1976 johnpye 669 /** create a term from the pool */
1977 aw0a 1 static struct relation_term *CreateRealTerm(double v, CONST dim_type *dim)
1978     {
1979     struct relation_term *term;
1980     term = POOL_ALLOCTERM;
1981     assert(term!=NULL);
1982     PTINIT(term);
1983     term->t = e_real;
1984     R_TERM(term) -> value = v;
1985     R_TERM(term) -> dimensions = dim;
1986     return term;
1987     }
1988    
1989 johnpye 669 /** create a term from the pool. Zero terms look like real, wild zeros */
1990 aw0a 1 static struct relation_term *CreateZeroTerm(void)
1991     {
1992     struct relation_term *term;
1993     term = POOL_ALLOCTERM;
1994     assert(term!=NULL);
1995     PTINIT(term);
1996     term->t = e_zero;
1997     R_TERM(term)->value = 0.0;
1998     R_TERM(term)->dimensions = WildDimension();
1999     return term;
2000     }
2001    
2002 johnpye 669 /** create a term from the pool */
2003 aw0a 1 static struct relation_term *CreateFuncTerm(CONST struct Func *f)
2004     {
2005     struct relation_term *term;
2006     term = POOL_ALLOCTERM;
2007     assert(term!=NULL);
2008     PTINIT(term);
2009     term->t = e_func;
2010     F_TERM(term) -> fptr = f;
2011     F_TERM(term) -> left = NULL;
2012     return term;
2013     }
2014    
2015 johnpye 669 /** create a term from the pool */
2016 aw0a 1 #ifdef THIS_IS_AN_UNUSED_FUNCTION
2017     static struct relation_term *CreateNaryTerm(CONST struct Func *f)
2018     {
2019     struct relation_term *term;
2020     term = POOL_ALLOCTERM;
2021     assert(term!=NULL);
2022     PTINIT(term);
2023     term->t = e_func;
2024 johnpye 709 N_TERM(term)->fptr = f;
2025     N_TERM(term)->args = NULL;
2026 aw0a 1 return term;
2027     }
2028     #endif /* THIS_IS_AN_UNUSED_FUNCTION */
2029    
2030    
2031 johnpye 669 /**
2032 johnpye 709 This function creates and *must* create the memory
2033     for the structure and for the union that the structure
2034 johnpye 908 points to.
2035 johnpye 709
2036     Too much code depends on the pre-existence of a properly initialized union.
2037    
2038     If copyunion is crs_NOUNION, the share ptr is init to NULL and user
2039     must set refcount,relop after the allocate a UNION or whatever.
2040     If copyunion is crs_NEWUNION, share ptr is allocated and configured.
2041     */
2042 aw0a 1 struct relation *CreateRelationStructure(enum Expr_enum relop,int copyunion)
2043     {
2044     struct relation *newrelation;
2045    
2046 johnpye 709 newrelation = ASC_NEW(struct relation);
2047 aw0a 1 assert(newrelation!=NULL);
2048 johnpye 731 /* CONSOLE_DEBUG("Created 'struct relation' at %p",newrelation); */
2049 aw0a 1
2050     newrelation->residual = DBL_MAX;
2051     newrelation->multiplier = DBL_MAX;
2052     newrelation->nominal = 1.0;
2053     newrelation->iscond = 0;
2054     newrelation->vars = NULL;
2055     newrelation->d =(dim_type *)WildDimension();
2056 johnpye 908 newrelation->externalData = NULL;
2057 aw0a 1
2058     if (copyunion) {
2059 johnpye 709 newrelation->share = ASC_NEW(union RelationUnion);
2060 aw0a 1 assert(newrelation->share!=NULL);
2061     RelationRefCount(newrelation) = 0;
2062     RelRelop(newrelation) = relop;
2063     #if TOKENDOMINANT
2064     RTOKEN(newrelation).lhs_term = NULL;
2065     RTOKEN(newrelation).rhs_term = NULL;
2066     RTOKEN(newrelation).lhs = NULL;
2067     RTOKEN(newrelation).rhs = NULL;
2068     RTOKEN(newrelation).lhs_len = 0;
2069     RTOKEN(newrelation).rhs_len = 0;
2070     RTOKEN(newrelation).btable = 0;
2071     RTOKEN(newrelation).bindex = 0;
2072     #else
2073     memset((char *)(newrelation->share),0,sizeof(union RelationUnion));
2074     #endif
2075     } else {
2076     newrelation->share = NULL;
2077     }
2078     return newrelation;
2079     }
2080    
2081 johnpye 669
2082     /*------------------------------------------------------------------------------
2083 johnpye 709 EXTERNAL CALL PROCESSING
2084     */
2085 johnpye 700
2086 johnpye 709 /** @file "relation.h"
2087 johnpye 908 @note
2088 johnpye 700 A special note on external relations
2089    
2090     External relations behave like relations but they also behave like
2091     procedures. As such when they are constructed and invoked they expect
2092     a particular ordering of their variables.
2093 johnpye 908
2094 johnpye 700 However there are some operations that can mess up (reduce) the number
2095     of incident variables on the incident varlist -- ATSing 2 variables in the
2096     *same* relation will do this. BUT we still need to maintain the number
2097     of variables in the call to the evaluation routine.
2098 johnpye 908
2099 johnpye 700 Consider the following example:
2100     An glassbox relation is constructed as: test1(x[46,2,8,9] ; 2);
2101     It *requires* 4 arguements, but its incident var count could be anything
2102 johnpye 908 from 1 <= n <= 4, depending on how many ATS are done.
2103     The ATS/alias will be done even before we have constructed the relation,
2104     so we just issue warnings.
2105 johnpye 669 */
2106 johnpye 700 struct relation *CreateBlackBoxRelation(struct Instance *relinst
2107     , struct Instance *subject
2108     , struct gl_list_t *inputs
2109 johnpye 908 , struct BlackBoxCache * common
2110     , unsigned long lhsIndex
2111     , CONST char *context
2112 johnpye 700 ){
2113 aw0a 1 struct relation *result;
2114 johnpye 908 struct gl_list_t *varlist;
2115     struct BlackBoxData *bbd;
2116 johnpye 89 struct Instance *var = NULL;
2117 johnpye 908 unsigned long *inputArgs;
2118 johnpye 1031 int32 inputsLen;
2119 johnpye 908 int argloc;
2120 johnpye 1031 unsigned long c,len,pos, lhsVarNumber;
2121 aw0a 1 unsigned long n_inputs;
2122 johnpye 908 CONST dim_type *d;
2123 aw0a 1
2124 jpye 1359 /* CONSOLE_DEBUG("CREATING BLACK BOX RELATION"); */
2125 johnpye 194
2126 aw0a 1 n_inputs = gl_length(inputs);
2127 johnpye 700 len = n_inputs + 1; /* an extra for the output variable. */
2128 aw0a 1
2129     /*
2130 johnpye 908 Create the BlackBox relation structure.
2131     output var always first in varlist.
2132 johnpye 700 */
2133 johnpye 1031 bbd = CreateBlackBoxData(common);
2134     inputsLen = BlackBoxCacheInputsLen(common);
2135     inputArgs = (unsigned long *)ascmalloc(sizeof(unsigned long) * inputsLen);
2136     lhsVarNumber = 1;
2137 johnpye 908 varlist = gl_create(len);
2138 johnpye 697
2139 johnpye 908 /* add the subject */
2140     gl_append_ptr(varlist,(VOIDPTR)subject); /* add the subject */
2141     AddRelation(subject,relinst);
2142    
2143     /* now loop, warning of merges and collecting varlist position
2144     of each arg into argloc.
2145     */
2146     argloc = 0;
2147     for (c=1; c<=n_inputs; c++) {
2148 aw0a 1 var = (struct Instance *)gl_fetch(inputs,c);
2149 johnpye 908 pos = gl_search(varlist,var,(CmpFunc)CmpP);
2150     switch (pos) {
2151     case 0:
2152     gl_append_ptr(varlist,(VOIDPTR)var);
2153 johnpye 725 AddRelation(var,relinst);
2154 johnpye 908 break;
2155     case 1:
2156 johnpye 1034 ERROR_REPORTER_HERE(ASC_USER_WARNING,"In external relation %s[%d],"
2157     " output %d and input %d are merged."
2158     " This will probably destroy your chances of achieving convergence,"
2159     " unless you are very careful."
2160     , context, (int)lhsIndex+1, (int)lhsIndex, (int)pos
2161     );
2162 johnpye 908 break;
2163     default:
2164 johnpye 1034 ERROR_REPORTER_HERE(ASC_USER_WARNING,"In external relation %s[%d],"
2165     " input %d and input %d are merged."
2166     " This will probably destroy your chances of achieving convergence,"
2167     " unless you are very careful."
2168     , context, (int)lhsIndex+1, argloc, (int)pos
2169     );
2170 johnpye 908 break;
2171 aw0a 1 }
2172 johnpye 908 if (pos) {
2173     inputArgs[argloc] = pos;
2174     } else {
2175     inputArgs[argloc] = gl_length(varlist);
2176     }
2177     argloc++;
2178 aw0a 1 }
2179    
2180 johnpye 709 /*
2181 johnpye 908 Now make the main relation structure and put it all
2182     together.
2183 johnpye 709 */
2184 johnpye 908 result = CreateRelationStructure(e_equal,crs_NEWUNION);
2185     RelationRefCount(result) = 1;
2186     result->externalData = bbd;
2187     result->vars = varlist;
2188 johnpye 1031 RBBOX(result).inputArgs = inputArgs;
2189     RBBOX(result).lhsindex = lhsIndex;
2190     RBBOX(result).lhsvar = lhsVarNumber;
2191 johnpye 908 d = RealAtomDims(subject);
2192     if (!IsWild(d)) {
2193     SetRelationDim(result,d);
2194 aw0a 1 }
2195 johnpye 908 return result;
2196     }
2197 aw0a 1
2198 johnpye 709
2199 aw0a 1 struct relation *CreateGlassBoxRelation(struct Instance *relinst,
2200     struct ExternalFunc *efunc,
2201     struct gl_list_t *varlist,
2202 johnpye 908 int gbindex,
2203 aw0a 1 enum Expr_enum relop)
2204     {
2205     struct relation *result;
2206     struct Instance *var;
2207     struct gl_list_t *newlist = NULL;
2208 johnpye 89 int *tmp = NULL, *args = NULL;
2209 aw0a 1 unsigned long len,c,pos;
2210    
2211     len = gl_length(varlist);
2212     /*
2213     * Make the variables aware that they are incident
2214     * in this relation instance. At the same time set up
2215     * the args list indexing.
2216     */
2217     if (len) {
2218 johnpye 709 tmp = args = ASC_NEW_ARRAY_CLEAR(int,len+1);
2219 aw0a 1 newlist = gl_create(len);
2220    
2221     for (c=1;c<=len;c++) {
2222     var = (struct Instance *)gl_fetch(varlist,c);
2223     pos = gl_search(newlist,var,(CmpFunc)CmpP);
2224     if (pos) {
2225     FPRINTF(ASCERR,"Incidence for external relation will be inaccurate\n");
2226     *tmp++ = (int)pos;
2227     }
2228     else{
2229     gl_append_ptr(newlist,(VOIDPTR)var);
2230     *tmp++ = (int)gl_length(newlist);
2231     AddRelation(var,relinst);
2232     }
2233     }
2234     }
2235 johnpye 700 *tmp = 0; /* terminate */
2236 aw0a 1
2237     /*
2238 johnpye 700 Create the relation data structure and append the
2239     varlist.
2240     */
2241 aw0a 1 result = CreateRelationStructure(relop,crs_NEWUNION);
2242     RelationRefCount(result) = 1;
2243     RGBOX(result).efunc = efunc;
2244     RGBOX(result).args = args;
2245     RGBOX(result).nargs = (int)len;
2246 johnpye 908 RGBOX(result).index = gbindex;
2247 aw0a 1 result->vars = newlist;
2248     return result;
2249     }
2250    
2251 johnpye 669 /*------------------------------------------------------------------------------
2252     TOKENRELATION PROCESSING AND GENERAL EXPR-TO-RELATION CHECK ROUTINES
2253     */
2254 aw0a 1
2255     static
2256     struct value_t CheckIntegerCoercion(struct value_t v)
2257     {
2258     if ((ValueKind(v)==real_value) && (RealValue(v)==0.0) &&
2259     IsWild(RealValueDimensions(v)) ){
2260     DestroyValue(&v);
2261     return CreateIntegerValue(0,1); /* assume this is a constant then */
2262     }
2263     else return v;
2264     }
2265    
2266     static
2267     int ProcessListRange(CONST struct Instance *ref,
2268     CONST struct Expr *low,
2269     CONST struct Expr *up,
2270     int *added,
2271     int i,
2272     enum relation_errors *err,
2273     enum find_errors *ferr)
2274     {
2275     struct value_t lower,upper;
2276     struct relation_term *term;
2277     long lv,uv;
2278     assert(GetEvaluationContext()==NULL);
2279     SetEvaluationContext(ref);
2280     lower = EvaluateExpr(low,NULL,InstanceEvaluateName);
2281     upper = EvaluateExpr(up,NULL,InstanceEvaluateName);
2282     SetEvaluationContext(NULL);
2283     lower = CheckIntegerCoercion(lower);
2284     upper = CheckIntegerCoercion(upper);
2285     if ((ValueKind(lower)==integer_value)&&(ValueKind(upper)==integer_value)){
2286     lv = IntegerValue(lower);
2287     uv = IntegerValue(upper);
2288     while(lv<=uv){
2289     term = CreateIntegerTerm(lv);
2290     AppendTermBuf(term);
2291     if ((*added)++) {
2292     switch(i){
2293     case SUM:
2294     term = CreateOpTerm(e_plus);
2295     break;
2296     case PROD:
2297     term = CreateOpTerm(e_times);
2298     break;
2299     }
2300     AppendTermBuf(term);
2301     }
2302     lv++;
2303     }
2304     return 0;
2305     }
2306     else{
2307     if(ValueKind(lower)==error_value) {
2308     FigureOutError(lower,err,ferr);
2309     return 1;
2310     }
2311     if(ValueKind(upper)==error_value){
2312     FigureOutError(upper,err,ferr);
2313     return 1;
2314     }
2315     *err = incorrect_structure;
2316     FPRINTF(ASCERR,"incorrect_structure in ProcessListRange\n");
2317     return 1;
2318     }
2319     }
2320    
2321     static
2322     CONST struct Expr *ExprContainsSuchThat(register CONST struct Expr *ex)
2323     {
2324     while(ex!=NULL){
2325     if (ExprType(ex)==e_st) return ex;
2326     ex = NextExpr(ex);
2327     }
2328     return ex;
2329     }
2330    
2331 johnpye 669 /**
2332 aw0a 1 * Here we give up if vars are not well defined.
2333     * At present e_var acceptable ARE:
2334     * REAL_ATOM_INSTANCE
2335     * Well defined Real and Integer constants.
2336     * Everything else is trash.
2337     * CreateTermFromInst() and CheckExpr() must have matching semantics.
2338     */
2339     static
2340     struct relation_term *CreateTermFromInst(struct Instance *inst,
2341     struct Instance *rel,
2342     enum relation_errors *err)
2343     {
2344     struct relation_term *term;
2345     switch(InstanceKind(inst)){
2346     case REAL_ATOM_INST:
2347     term = CreateVarTerm(inst);
2348     AddRelation(inst,rel);
2349     return term;
2350     case REAL_CONSTANT_INST:
2351     if ( AtomAssigned(inst) && !IsWild(RealAtomDims(inst)) ){
2352     term = CreateRealTerm(RealAtomValue(inst),RealAtomDims(inst));
2353     return term;
2354     }
2355     else{
2356     if ( IsWild(RealAtomDims(inst)) && AtomAssigned(inst) ) {
2357     *err = real_value_wild;
2358     } else {
2359     *err = real_value_undefined;
2360     }
2361     return NULL;
2362     }
2363     case INTEGER_CONSTANT_INST:
2364     if (AtomAssigned(inst)){
2365     term = CreateIntegerTerm(GetIntegerAtomValue(inst));
2366     return term;
2367     }
2368     else{
2369     *err = integer_value_undefined;
2370     return NULL;
2371     }
2372     case REAL_INST:
2373     *err = incorrect_real_inst_type;
2374     return NULL;
2375     case INTEGER_ATOM_INST:
2376     case INTEGER_INST:
2377     *err = incorrect_integer_inst_type;
2378     return NULL;
2379     case SYMBOL_ATOM_INST:
2380     case SYMBOL_CONSTANT_INST:
2381     case SYMBOL_INST:
2382     *err = incorrect_symbol_inst_type;
2383     return NULL;
2384     case BOOLEAN_ATOM_INST:
2385     case BOOLEAN_CONSTANT_INST:
2386     case BOOLEAN_INST:
2387     *err = incorrect_boolean_inst_type;
2388     return NULL;
2389     default:
2390     *err = incorrect_inst_type;
2391     return NULL;
2392     }
2393     }
2394    
2395     /* forward declaration */
2396     static int AppendList( CONST struct Instance *,
2397     struct Instance *,
2398     CONST struct Set *,
2399     int ,
2400     enum relation_errors *,
2401     enum find_errors *);
2402    
2403 johnpye 669 /**
2404     @todo document this
2405 johnpye 908
2406 johnpye 669 Convert a part of an expression into part of a relation (in postfix)?
2407     */
2408 aw0a 1 static
2409     int ConvertSubExpr(CONST struct Expr *ptr,
2410     CONST struct Expr *stop,
2411     CONST struct Instance *ref,
2412     struct Instance *rel,
2413     int *added,
2414     int i,
2415     enum relation_errors *err,
2416     enum find_errors *ferr)
2417     {
2418 johnpye 89 struct relation_term *term = NULL;
2419 aw0a 1 struct gl_list_t *instances;
2420     unsigned c,len;
2421     struct Instance *inst;
2422     struct value_t svalue,cvalue;
2423     int my_added=0;
2424     symchar *str;
2425     CONST struct for_var_t *fvp; /* for var pointer */
2426     while (ptr!=stop){
2427     switch(ExprType(ptr)){
2428     case e_plus:
2429     case e_minus:
2430     case e_times:
2431     case e_divide:
2432     case e_power:
2433     case e_ipower:
2434     case e_uminus:
2435     term = CreateOpTerm(ExprType(ptr));
2436     my_added++;
2437     AppendTermBuf(term);
2438     break;
2439     case e_var:
2440     str = SimpleNameIdPtr(ExprName(ptr));
2441     if (str&&TempExists(str)){
2442 johnpye 89 cvalue = TempValue(str);
2443     switch(ValueKind(cvalue)){
2444     case integer_value:
2445     term = CreateIntegerTerm(IntegerValue(cvalue));
2446     my_added++;
2447     AppendTermBuf(term);
2448     break;
2449     default:
2450     FPRINTF(ASCERR,"Non-integer temporary variable used in expression.\n");
2451     *err = incorrect_inst_type;
2452     term = NULL;
2453     return 1;
2454     }
2455     }else if (GetEvaluationForTable() != NULL && str !=NULL &&
2456 aw0a 1 (fvp=FindForVar(GetEvaluationForTable(),str)) !=NULL ){
2457 johnpye 89 if (GetForKind(fvp)==f_integer){
2458     term = CreateIntegerTerm(GetForInteger(fvp));
2459     my_added++;
2460     AppendTermBuf(term);
2461     }
2462     else{
2463     FPRINTF(ASCERR,
2464 aw0a 1 "Non-integer FOR variable used in expression.\n");
2465 johnpye 89 *err = incorrect_inst_type;
2466     return 1;
2467     }
2468 aw0a 1 }
2469     else{
2470 johnpye 669 instances = FindInstances(ref,ExprName(ptr),ferr);
2471     if (instances!=NULL){
2472     if (NextExpr(ptr)==stop){ /* possibly multiple instances */
2473     len = gl_length(instances);
2474     for(c=1;c<=len;c++){
2475     inst = (struct Instance *)gl_fetch(instances,c);
2476     if ((term=CreateTermFromInst(inst,rel,err))!=NULL){
2477     AppendTermBuf(term);
2478     if (my_added++){
2479     switch(i){
2480     case SUM:
2481     term = CreateOpTerm(e_plus);
2482     break;
2483     case PROD:
2484     term = CreateOpTerm(e_times);
2485     break;
2486     }
2487     AppendTermBuf(term);
2488     }
2489     }
2490     else{
2491     gl_destroy(instances);
2492     return 1;
2493     }
2494     }
2495     gl_destroy(instances);
2496     }
2497     else{ /* single instance */
2498     if (gl_length(instances)==1){
2499     inst = (struct Instance *)gl_fetch(instances,1);
2500     gl_destroy(instances);
2501     if ((term=CreateTermFromInst(inst,rel,err))!=NULL){
2502     my_added++;
2503     AppendTermBuf(term);
2504     }
2505     else
2506     return 1;
2507     }
2508     else{
2509     gl_destroy(instances);
2510     *err = incorrect_structure;
2511 aw0a 1 FPRINTF(ASCERR,"incorrect_structure in ConvertSubExpr 1\n");
2512 johnpye 669 return 1;
2513     }
2514     }
2515     } else{
2516     *err = find_error;
2517     return 1;
2518     }
2519 aw0a 1 }
2520     break;
2521     case e_int:
2522     term = CreateIntegerTerm(ExprIValue(ptr));
2523     my_added++;
2524     AppendTermBuf(term);
2525     break;
2526     case e_zero:
2527     /* this should never happen here */
2528     term = CreateZeroTerm();
2529     my_added++;
2530     AppendTermBuf(term);
2531     break;
2532     case e_real:
2533     term = CreateRealTerm(ExprRValue(ptr),ExprRDimensions(ptr));
2534     my_added++;
2535     AppendTermBuf(term);
2536     break;
2537     case e_card:
2538     assert(GetEvaluationContext() == NULL);
2539     SetEvaluationContext(ref);
2540     svalue = EvaluateSet(ExprBuiltinSet(ptr),InstanceEvaluateName);
2541     SetEvaluationContext(NULL);
2542     cvalue = CardValues(svalue);
2543     DestroyValue(&svalue);
2544     switch(ValueKind(cvalue)){
2545     case integer_value:
2546 johnpye 669 term = CreateIntegerTerm(IntegerValue(cvalue));
2547     my_added++;
2548     AppendTermBuf(term);
2549     break;
2550 aw0a 1 case error_value:
2551 johnpye 669 FigureOutError(cvalue,err,ferr);