/[ascend]/trunk/base/generic/compiler/exprsym.c
ViewVC logotype

Contents of /trunk/base/generic/compiler/exprsym.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 423 - (show annotations) (download) (as text)
Mon Apr 3 23:15:49 2006 UTC (18 years, 2 months ago) by ben.allan
File MIME type: text/x-csrc
File size: 31303 byte(s)
Fixed instatiation of blackbox (and probably glassbox)
relations, for poorly tested definitions of 'fix'.
Blackbox evaluation is still broken.

There seems to be some chaos around win32_lean_and_mean
in ascConfig.h

updated reconfig. setenv TKROOT /where/is/tk8.3dir before 
running and all is good.
1 /*
2 * Symbolic Expression Manipulation
3 * by Kirk Abbott
4 * Created: Novermber 21, 1994
5 * Version: $Revision: 1.9 $
6 * Version control file: $RCSfile: exprsym.c,v $
7 * Date last modified: $Date: 1997/09/08 18:07:36 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the ASCEND compiler.
11 *
12 * Copyright (C) 1994,1995 Kirk Andre Abbott.
13 *
14 * The Ascend Language Interpreter is free software; you can redistribute
15 * it and/or modify it under the terms of the GNU General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * The Ascend Language Interpreter is distributed in hope that it will be
20 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with the program; if not, write to the Free Software Foundation,
26 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
27 * COPYING. COPYING is found in ../compiler.
28 *
29 */
30
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <math.h> /* !!! dont ever forget this and -lm */
34 #include <utilities/ascConfig.h>
35 #include "compiler.h"
36 #include <utilities/ascMalloc.h>
37 #include <general/list.h>
38 #include "fractions.h"
39 #include "dimen.h"
40 #include "functype.h"
41 #include "func.h"
42 #include "types.h"
43 #include "instance_enum.h"
44 #include "freestore.h"
45 #include "relation_type.h"
46 #include "find.h"
47 #include "relation.h"
48 #include "relation_util.h"
49 #include "exprsym.h"
50
51 #ifndef NO_FREE_STORE
52 /* free store alloc */
53 #define TMALLOC (Term *)GetMem()
54 /* independent allocation */
55 #define UTMALLOC TERM_ALLOC
56 #else
57 /* no free store --> both are independent allocation */
58 #define TMALLOC TERM_ALLOC
59 #define UTMALLOC TERM_ALLOC
60 #endif
61 #define RMALLOC (RelationINF *)ascmalloc((unsigned)sizeof(RelationINF))
62
63 /*
64 * Toggle the next 2 for some memory usage statistics.
65 */
66 #define FREESTORE_STATS 1
67 #undef FREESTORE_STATS
68 static int CmpP(CONST char *c1, CONST char *c2)
69 {
70 if (c1 > c2) return 1;
71 if (c1 < c2) return -1;
72 return 0;
73 }
74
75 /*
76 * Very popular constants. We have not put them in the free store
77 * as it takes too much fancy footwork not to overwrite them each
78 * time that we (re)init the free store.
79 */
80 /* note that Relation real is the size dominant relation_term
81 subtype, so these are ok.
82 */
83 static struct RelationReal MinusOne;
84 static struct RelationReal Zero;
85 static struct RelationReal One;
86 static struct RelationReal Two;
87 static struct RelationReal Three;
88 static int constants_inited = 0;
89
90 /*********************************************************************\
91 Terms, and routines to manipulate them.
92
93 At the moment the orginal relations use malloc for getting
94 memory, rather than making use of a freestore. Consequently
95 we have *very* similar code for making the Terms. One day the
96 the original relations will use a freestore as well in which
97 case, the term creation routines may be unified. In the mean time
98 all the Create routines have a "D" (for derivative) inserted in their
99 names and have been made static.
100
101 Actually, however, the free stores, or pool stores, must be separate
102 to maintain memory integrity, so the create functions in
103 relation.c are private. These are free to become global if that
104 is desired. Probably not desired however...
105 \*********************************************************************/
106
107 static Term *CreateDTermBinary(Term *left, enum Expr_enum kind, Term *right)
108 {
109 Term *result;
110 result = TMALLOC;
111 result->t = kind; /* set the kind of node */
112 B_TERM(result)->left = left;
113 B_TERM(result)->right = right;
114 return result;
115 }
116
117 #ifdef THIS_IS_AN_UNUSED_FUNCTION
118 static Term *CreateDTermUnary(Term *left, enum Expr_enum kind)
119 {
120 Term *result;
121 result = TMALLOC;
122 result->t = kind; /* set the kind of node */
123 U_TERM(result)->left = left;
124 return result;
125 }
126 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
127
128
129 #ifdef THIS_IS_AN_UNUSED_FUNCTION
130 static Term *CreateDTermVar(void)
131 {
132 Term *result;
133 result = TMALLOC;
134 result->t = e_var;
135 V_TERM(result)->varnum = 0; /* not given a number yet */
136 return result;
137 }
138 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
139
140
141 static Term *CreateDTermInteger(long int i)
142 {
143 Term *result;
144 result = TMALLOC;
145 result->t = e_int;
146 I_TERM(result)->ivalue = i;
147 return result;
148 }
149
150 static Term *CreateDTermReal(double d, dim_type *dim)
151 {
152 Term *result;
153 result = TMALLOC;
154 result->t = e_real;
155 R_TERM(result)->value = d;
156 R_TERM(result)->dimensions = dim;
157 return result;
158 }
159
160 static Term *CreateDTermFunc(Term *left,CONST Func *fptr)
161 {
162 Term *result;
163 result = TMALLOC;
164 result->t = e_func;
165 F_TERM(result)->left = left;
166 F_TERM(result)->fptr = (Func *)fptr;
167 return result;
168 }
169
170
171 /*
172 * Clone and Copy degenerate into the same
173 * function whenever, the freestore is not in
174 * use. Copy must be used whenever a caller
175 * requires his own *separate* copy of a relation
176 * and its associated terms.
177 */
178 /* pull a copy from the free store */
179 static Term *CloneDTerm(Term *term)
180 {
181 Term *result;
182 result = TMALLOC;
183 memcpy((VOIDPTR)result,(VOIDPTR)term,sizeof(union RelationTermUnion));
184 return result;
185 }
186
187 /* create an independent copy */
188 static Term *CopyDTerm(Term *term)
189 {
190 Term *result;
191
192 assert(term!=NULL);
193 result = UTMALLOC; /* independent copy */
194 memcpy((VOIDPTR)result,(VOIDPTR)term,sizeof(union RelationTermUnion));
195 return result;
196 }
197
198 #ifdef THIS_IS_AN_UNUSED_FUNCTION
199 static void DestroyDTerm(Term *term)
200 {
201 FreeMem(UNION_TERM(term));
202 }
203 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
204
205 /* init static constant real terms we use a lot. must remember to cast
206 * them appropriately elsewhere
207 */
208 static void InitConstantTerms(void)
209 {
210 if (!constants_inited) {
211 MinusOne.t = e_real;
212 MinusOne.value = -1.0;
213 MinusOne.dimensions = NULL; /* should probably make wild */
214
215 Zero.t = e_real;
216 Zero.value = 0.0;
217 Zero.dimensions = NULL;
218
219 One.t = e_real;
220 One.value = 1.0;
221 One.dimensions = NULL;
222
223 Two.t = e_real;
224 Two.value = 2.0;
225 Two.dimensions = NULL;
226
227 Three.t = e_real;
228 Three.value = 3.0;
229 Three.dimensions = NULL;
230 constants_inited = 1;
231 }
232 }
233
234 /* what the heck is this here for? */
235 static
236 RelationINF *CreateRelInfix(Term *lhs, enum Expr_enum op, Term *rhs)
237 {
238 RelationINF *result;
239 result = (RelationINF *)CreateRelationStructure(op,crs_NEWUNION);
240 result->share->s.relop = op;
241 RelationRefCount(result) = 1;
242 RTOKEN(result).lhs = NULL; /* this is for postfix scanning */
243 RTOKEN(result).rhs = NULL;
244 RTOKEN(result).lhs_term = lhs;
245 RTOKEN(result).rhs_term = rhs;
246 result->vars = NULL;
247 result->d = NULL;
248 result->residual = 0.0;
249 result->multiplier = 0.0;
250 return result;
251 }
252
253 /*
254 * This function will run through a tree of terms
255 * and fix the varlist associated with the tree.
256 * It *assumes* that the current tree has a lower
257 * # of inident variables that the old list.
258 */
259 static struct gl_list_t *g_new_var_list = NULL;
260
261 static
262 void CollectVars(Term *term, RelationINF *rel)
263 {
264 unsigned long oldvarnum, pos;
265 CONST struct Instance *inst;
266
267 assert(term!=NULL);
268 switch(RelationTermType(term)) {
269 case e_var:
270 oldvarnum = TermVarNumber(term);
271 inst = RelationVariable(rel,oldvarnum);
272 pos = gl_search(g_new_var_list,inst,(CmpFunc)CmpP);
273 if (pos)
274 V_TERM(term)->varnum = pos;
275 else{
276 gl_append_ptr(g_new_var_list,(VOIDPTR)inst);
277 V_TERM(term)->varnum = gl_length(g_new_var_list);
278 }
279 break;
280 default:
281 break;
282 }
283 }
284
285 static
286 struct gl_list_t *RedoVarList(RelationINF *rel,
287 struct gl_list_t *oldvars)
288 {
289 unsigned long oldlen;
290 struct gl_list_t *tmp;
291 oldlen = (oldvars!=NULL) ? gl_length(oldvars) : 0;
292 if (oldlen) {
293 g_new_var_list = gl_create(oldlen); /* be conservative ! */
294 DoInOrderVisit(Infix_LhsSide(rel),rel,CollectVars);
295 DoInOrderVisit(Infix_RhsSide(rel),rel,CollectVars);
296 tmp = g_new_var_list;
297 g_new_var_list = NULL;
298 return tmp;
299 }
300 return NULL;
301 }
302
303 /*
304 * Return a separate verbatim copy of a term tree.
305 * i.e, no simplification, etc..
306 */
307
308 static
309 Term *CopyDTermTree(Term *term)
310 {
311 Term *left,*right,*result;
312 if (!term) return NULL;
313 switch(term->t) {
314 case e_var: case e_int: case e_real: case e_zero:
315 return CopyDTerm(term);
316 case e_plus: case e_minus: case e_times:
317 case e_divide: case e_power: case e_ipower:
318 left = CopyDTermTree(TermBinLeft(term));
319 right = CopyDTermTree(TermBinRight(term));
320 result = CopyDTerm(term);
321 B_TERM(result)->left = left; B_TERM(result)->right = right;
322 return result;
323 case e_uminus:
324 left = CopyDTermTree(TermUniLeft(term));
325 result = CopyDTerm(term);
326 U_TERM(result)->left = left;
327 return result;
328 case e_func:
329 left = CopyDTermTree(TermFuncLeft(term));
330 result = CopyDTerm(term);
331 F_TERM(result)->left = left;
332 return result;
333 default:
334 FPRINTF(ASCERR,"Unknown term type in expression\n");
335 return NULL;
336 }
337 }
338
339
340 /*********************************************************************\
341 Support routines for symbolic manipulation
342 \*********************************************************************/
343
344 /*
345 * The basic principle here is that all the Make* functions
346 * must *never* modify what has been given to them. They
347 * must create memory and return it, if necessary.
348 */
349 static Term *MakeAdd(Term *left, Term *right)
350 {
351 Term *result = NULL;
352
353 if (left == K_TERM(&Zero)) return right; /* 0 + u = u */
354 if (right == K_TERM(&Zero)) return left; /* u + 0 = u */
355
356 if (left->t == e_real && right->t == e_real)
357 return CreateDTermReal(R_TERM(left)->value + R_TERM(right)->value, NULL);
358
359 if (left->t == e_real && right->t == e_int)
360 return CreateDTermReal(R_TERM(left)->value
361 + (double)I_TERM(right)->ivalue, NULL);
362
363 if (left->t == e_int && right->t == e_real)
364 return CreateDTermReal((double)I_TERM(left)->ivalue
365 + R_TERM(right)->value, NULL);
366
367 if (left->t == e_int && right->t == e_int)
368 return CreateDTermInteger(I_TERM(left)->ivalue + I_TERM(right)->ivalue);
369
370 result = TMALLOC;
371 if (right->t == e_uminus) { /* a + (-b) --> a - b */
372 result->t = e_minus;
373 B_TERM(result)->left = left;
374 B_TERM(result)->right = TermUniLeft(right);
375 return result;
376 }
377 if (left->t == e_uminus) { /* (-a) + b --> b - a */
378 result->t = e_minus;
379 B_TERM(result)->left = right;
380 B_TERM(result)->right = TermUniLeft(left);
381 }
382
383 result->t = e_plus;
384 B_TERM(result)->left = left;
385 B_TERM(result)->right = right;
386 return result;
387 }
388
389 static Term *MakeSubtract(Term *left, Term *right)
390 {
391 Term *result=NULL;
392
393 if (right==K_TERM(&Zero)) return left; /* u - 0 = u */
394
395 if (left->t == e_real && right->t == e_real)
396 return CreateDTermReal(R_TERM(left)->value - R_TERM(right)->value, NULL);
397
398 if (left->t == e_real && right->t == e_int)
399 return CreateDTermReal(R_TERM(left)->value
400 - (double)I_TERM(right)->ivalue, NULL);
401
402 if (left->t == e_int && right->t == e_real)
403 return CreateDTermReal((double)I_TERM(left)->ivalue
404 - R_TERM(right)->value, NULL);
405
406 if (left->t == e_int && right->t == e_int)
407 return CreateDTermInteger(I_TERM(left)->ivalue - I_TERM(right)->ivalue);
408
409
410 result = TMALLOC;
411 if (right->t == e_uminus) {
412 result->t = e_plus; /* a - (-b) --> a + b */
413 B_TERM(result)->left = left;
414 B_TERM(result)->right = TermUniLeft(right);
415 return result;
416 }
417
418 result->t = e_minus;
419 B_TERM(result)->left = left;
420 B_TERM(result)->right = right;
421 return result;
422 }
423
424 static Term *MakeMultiply(Term *left, Term *right)
425 {
426 Term *result;
427
428 if (left==K_TERM(&Zero) || right == K_TERM(&Zero))
429 return K_TERM(&Zero);
430 if (left==K_TERM(&One) && right==K_TERM(&One)) return K_TERM(&One);
431 if (left==K_TERM(&One) && right!=K_TERM(&One)) return right;
432 if (left!=K_TERM(&One) && right==K_TERM(&One)) return left;
433
434 /* the following assumes no &One in the expression */
435
436 if (left->t == e_real && right->t == e_real)
437 return CreateDTermReal(R_TERM(left)->value * R_TERM(right)->value, NULL);
438
439 if (left->t == e_real && right->t == e_int)
440 return CreateDTermReal(R_TERM(left)->value * (double)I_TERM(right)->ivalue,
441 NULL);
442
443 if (left->t == e_int && right->t == e_real)
444 return CreateDTermReal((double)I_TERM(left)->ivalue * R_TERM(right)->value,
445 NULL);
446
447 if (left->t == e_int && right->t == e_int)
448 return CreateDTermInteger(I_TERM(left)->ivalue * I_TERM(right)->ivalue);
449
450 /* not a simple constant */
451
452 result = TMALLOC;
453 result->t = e_times;
454 B_TERM(result)->left = left;
455 B_TERM(result)->right = right;
456 return result;
457 }
458
459
460 static Term *MakeDivide(Term *left, Term *right)
461 {
462 Term *result;
463
464 if (left==K_TERM(&Zero)) return K_TERM(&Zero); /* 0 / u = 0 */
465 if (right==K_TERM(&Zero)) {
466 FPRINTF(ASCERR,"Potential for floating point dvision in MakeDivide\n");
467 /* continue processing. */
468 }
469 if (left==K_TERM(&One) && right==K_TERM(&One))
470 return K_TERM(&One);
471 if (left!=K_TERM(&One) && right==K_TERM(&One))
472 return left;
473
474 if (left->t == e_real && right->t == e_real)
475 return CreateDTermReal(R_TERM(left)->value / R_TERM(right)->value, NULL);
476
477 if (left->t == e_real && right->t == e_int)
478 return CreateDTermReal(R_TERM(left)->value / (double)I_TERM(right)->ivalue,
479 NULL);
480
481 if (left->t == e_int && right->t == e_real)
482 return CreateDTermReal((double)I_TERM(left)->ivalue / R_TERM(right)->value,
483 NULL);
484
485 result = TMALLOC;
486 result->t = e_divide;
487 B_TERM(result)->left = left;
488 B_TERM(result)->right = right;
489 return result;
490 }
491
492 static Term *MakePower(Term *left, Term *right)
493 {
494 Term *result;
495 enum Expr_enum ltype, rtype;
496 double lvalue,rvalue;
497 ltype = left->t; rtype = right->t;
498
499 if (left==K_TERM(&One)) /* 1 ^ ? = 1 */
500 return K_TERM(&One);
501 if (right==K_TERM(&One)) /* ? ^ 1 = ? */
502 return left;
503 if (rtype == e_int && I_TERM(right)->ivalue == 0) /* ? ^ 0 = 1 */
504 return K_TERM(&One);
505
506 if (ltype == e_real && rtype== e_real)
507 return CreateDTermReal(pow(R_TERM(left)->value,R_TERM(right)->value),NULL);
508
509 if (ltype == e_real && rtype == e_int) {
510 rvalue = (double)I_TERM(right)->ivalue;
511 return CreateDTermReal(pow(R_TERM(left)->value,rvalue), NULL);
512 }
513 if (ltype == e_int && rtype == e_real) {
514 lvalue = (double)R_TERM(left)->value;
515 return CreateDTermReal(pow(lvalue,R_TERM(right)->value), NULL);
516 }
517 if (ltype == e_int && rtype == e_int) {
518 lvalue = (double)I_TERM(left)->ivalue;
519 rvalue = (double)I_TERM(right)->ivalue;
520 return CreateDTermReal(pow(lvalue,rvalue),NULL);
521 }
522 result = TMALLOC;
523 result->t = e_power;
524 B_TERM(result)->left = left;
525 B_TERM(result)->right = right;
526 return result;
527 }
528
529 static Term *MakeUMinus(Term *left)
530 {
531 Term *result;
532 result = TMALLOC;
533 result->t = e_uminus;
534 U_TERM(result)->left = left;
535 return result;
536 }
537
538 static Term *MakeNegation(Term *term)
539 {
540 Term *result;
541 if (term==K_TERM(&Zero))
542 return term;
543 if (term==K_TERM(&One))
544 return K_TERM(&MinusOne);
545 if(term==K_TERM(&MinusOne))
546 return K_TERM(&One);
547 /*
548 * We cannot just compute the (-1.0 * value) as we might
549 * be passed in one of the special constants, and we *must*
550 * not corrupt those constants.
551 * We dont want to search over all possible special
552 * constants, so it is easier to create a new node.
553 */
554 switch(term->t) {
555 case e_real:
556 result = TMALLOC;
557 result->t = e_real;
558 R_TERM(result)->value = -1.0*R_TERM(term)->value;
559 return result;
560 case e_int:
561 result = TMALLOC;
562 result->t = e_int;
563 I_TERM(result)->ivalue = -1*I_TERM(term)->ivalue;
564 return result;
565 case e_uminus:
566 result = TermUniLeft(term);
567 return result;
568 default:
569 result = TMALLOC;
570 result->t = e_uminus;
571 U_TERM(result)->left = term;
572 return result;
573 }
574 }
575
576 #ifdef THIS_IS_AN_UNUSED_FUNCTION
577 static void DoBreakPoint(void)
578 {
579 return;
580 }
581 #endif /* THIS_IS_AN_UNUSED_FUNCTION */
582
583 static Term *MakeFunc(enum Func_enum id, Term *term)
584 {
585 Term *result;
586 CONST struct Func *fptr;
587 double value;
588
589 fptr = LookupFuncById(id);
590 assert(fptr!=NULL);
591 if (term->t == e_real) {/* evaluate the term */
592 value = FuncEval(fptr,R_TERM(term)->value);
593 result = CreateDTermReal(value,NULL);
594 return result;
595 }
596
597 if (term->t == e_int) {/* evaluate the term */
598 value = FuncEval(fptr,(double)I_TERM(term)->ivalue);
599 result = CreateDTermReal(value,NULL);
600 return result;
601 }
602
603 if (term->t == e_uminus) {
604 switch(id) {
605 case F_SIN: case F_SINH: /* these are the odd functions */
606 case F_TAN: case F_TANH: /* sin (-x) = - sin (x) */
607 case F_ARCSIN: case F_ARCSINH:
608 case F_ARCTAN: case F_ARCTANH:
609 #ifdef HAVE_ERF
610 case F_ERF:
611 #endif /* HAVE_ERF */
612 case F_CUBE:
613 case F_CBRT:
614 result = CreateDTermFunc(TermUniLeft(term),fptr);
615 result = MakeUMinus(result);
616 return result;
617 case F_COS: case F_COSH: /* these are the even functions */
618 case F_ARCCOS: case F_SQR: /* cos (-x) = cos(x) */
619 result = CreateDTermFunc(TermUniLeft(term),fptr);
620 return result;
621 default:
622 break;
623 }
624 }
625 result = CreateDTermFunc(term,fptr);
626 return result;
627 }
628
629 Term *TermSimplify(Term *term)
630 {
631 Term *left,*right;
632 if (term) {
633 switch(term->t) {
634 case e_var:
635 case e_int:
636 case e_real:
637 return term;
638 case e_plus:
639 left = TermSimplify(TermBinRight(term));
640 right = TermSimplify(TermBinRight(term));
641 return MakeAdd(left,right);
642 case e_minus:
643 left = TermSimplify(TermBinRight(term));
644 right = TermSimplify(TermBinRight(term));
645 return MakeSubtract(left,right);
646 case e_times:
647 left = TermSimplify(TermBinRight(term));
648 right = TermSimplify(TermBinRight(term));
649 return MakeMultiply(left,right);
650 case e_divide:
651 left = TermSimplify(TermBinRight(term));
652 right = TermSimplify(TermBinRight(term));
653 return MakeDivide(left,right);
654 case e_power:
655 case e_ipower:
656 left = TermSimplify(TermBinRight(term));
657 right = TermSimplify(TermBinRight(term));
658 return MakePower(left,right);
659 case e_uminus:
660 left = TermSimplify(TermUniLeft(term));
661 return MakeUMinus(left);
662 case e_func:
663 left = TermSimplify(TermFuncLeft(term));
664 return MakeFunc(FuncId(TermFunc(term)),left);
665 default:
666 FPRINTF(ASCERR,"Unknown term type in expression\n");
667 return NULL;
668 }
669 }
670 return NULL;
671 }
672
673
674
675 /*********************************************************************\
676 Symbolic derivatives.
677 \*********************************************************************/
678
679 /*
680 * Any potentially destructive operations should be preceeded
681 * by a copying of the terms. This is necessary so that we don't
682 * trash the original expression tree that is handed to us.
683 * Still not sure whether I should have the filter operate on Terms
684 * or on instances. If it operated on terms then there would be the
685 * oppurtunity to handle reductions of e_int to e_real etc. For the
686 * moment it does nothing !!
687 */
688
689 Term *Derivative(Term *term, unsigned long wrt,
690 int (*filter)(struct Instance *))
691 {
692 Term *du_dx, *dv_dx;
693
694 switch(term->t) {
695 case e_var: /* du/dx = du_dx */
696 if (wrt==TermVarNumber(term))
697 return K_TERM(&One);
698 return K_TERM(&Zero);
699
700 case e_int:
701 case e_real:
702 return K_TERM(&Zero);
703
704 case e_plus: /* u + v */
705 du_dx = Derivative(TermBinRight(term),wrt,filter);
706 dv_dx = Derivative(TermBinRight(term),wrt,filter);
707
708 if (du_dx==K_TERM(&Zero) && dv_dx!=K_TERM(&Zero))
709 return dv_dx;
710 if (dv_dx==K_TERM(&Zero) && du_dx!=K_TERM(&Zero))
711 return du_dx;
712 if (du_dx==K_TERM(&Zero) && dv_dx==K_TERM(&Zero))
713 return K_TERM(&Zero);
714 return MakeAdd(du_dx,dv_dx);
715
716 case e_minus: /* u - v */
717 du_dx = Derivative(TermBinRight(term),wrt,filter);
718 dv_dx = Derivative(TermBinRight(term),wrt,filter);
719
720 if (du_dx==K_TERM(&Zero) && dv_dx!=K_TERM(&Zero))
721 return MakeNegation(dv_dx);
722 if (dv_dx==K_TERM(&Zero) && du_dx!=K_TERM(&Zero))
723 return du_dx;
724 if (du_dx==K_TERM(&Zero) && dv_dx==K_TERM(&Zero))
725 return K_TERM(&Zero);
726
727 return MakeSubtract(du_dx,dv_dx);
728
729 case e_times: { /* d(u * v)/dx --> (v * du_dx) + (u * dv_dx) */
730 Term *left, *right;
731
732 left = CloneDTerm(TermBinRight(term)); /* copy first !! */
733 right = CloneDTerm(TermBinRight(term));
734 du_dx = Derivative(left,wrt,filter);
735 dv_dx = Derivative(right,wrt,filter);
736
737 if (du_dx==K_TERM(&Zero) && dv_dx==K_TERM(&Zero)) /* --> 0 */
738 return K_TERM(&Zero);
739
740 if (du_dx==K_TERM(&Zero) && dv_dx!=K_TERM(&Zero)) /* --> u * dv_dx */
741 return MakeMultiply(left,dv_dx);
742
743 if (dv_dx==K_TERM(&Zero) && du_dx!=K_TERM(&Zero)) /* --> v * du_dx */
744 return MakeMultiply(right,du_dx);
745
746 return MakeAdd(MakeMultiply(left,dv_dx),
747 MakeMultiply(right,du_dx));
748 }
749 case e_divide: {/* d(u/v)/dx = (v*du_dx - u*dv_dx)/ v^2 */
750 /* = (du_dx/v - u/v^2*dv_dx) */
751 Term *left, *right;
752
753 left = CloneDTerm(TermBinRight(term));
754 right = CloneDTerm(TermBinRight(term));
755 du_dx = Derivative(left,wrt,filter);
756 dv_dx = Derivative(right,wrt,filter);
757
758 if (du_dx==K_TERM(&Zero) && dv_dx==K_TERM(&Zero)) /* --> 0 */
759 return K_TERM(&Zero);
760
761 if (du_dx==K_TERM(&Zero) && dv_dx!=K_TERM(&Zero)) {
762 /* --> - (u*dv_dx)/(v^2) */
763 return MakeNegation(MakeDivide(MakeMultiply(left,dv_dx),
764 MakeFunc(F_SQR,right)));
765 }
766
767 if (dv_dx==K_TERM(&Zero) && du_dx!=K_TERM(&Zero)) /* --> du_dx / v */
768 return MakeDivide(du_dx,right);
769
770 return MakeSubtract(MakeDivide(du_dx,right),
771 MakeDivide(MakeMultiply(left,dv_dx),
772 MakeFunc(F_SQR,right)));
773 }
774
775 case e_power:
776 case e_ipower: {/* d ( u ^ v)/dx = u^v * (du_dx * v / u + dv_dx*ln(u)) */
777 Term *left, *right;
778 Term *e0, *e1;
779 long ivalue;
780 int sival;
781 double rvalue;
782
783 left = CloneDTerm(TermBinRight(term));
784 right = CloneDTerm(TermBinRight(term));
785
786 /* 0 ^0 is undefined, so return Zero */
787 if (left==K_TERM(&Zero) && right==K_TERM(&Zero)) /* d (0^0)/dx = 0 */
788 return K_TERM(&Zero);
789 if (right==K_TERM(&Zero)) /* d (u^0)/dx = d(1)/dx = 0 */
790 return K_TERM(&Zero);
791 if (right==K_TERM(&One)) /* d (u^1)/dx = du_dx */
792 return Derivative(left,wrt,filter);
793 if (left==K_TERM(&One)) /* d(1^u)/dx = d(1)/dx = 0 */
794 return K_TERM(&Zero);
795
796 if (right->t==e_int) {
797 ivalue = I_TERM(right)->ivalue;
798 sival = (int)ivalue;
799 switch(sival) {
800 case 0: /* d (u^0)/dx : as above but comparison by value */
801 return K_TERM(&Zero);
802 case 1: /* d (u^1)/dx : as above but comparison by value */
803 return Derivative(left,wrt,filter);
804 case 2: /* d/dx (u ^ 2) = 2 * u * du_dx */
805 du_dx = Derivative(left,wrt,filter);
806 e0 = MakeMultiply(K_TERM(&Two),MakeMultiply(left,du_dx));
807 return e0;
808 default: /* d/dx (u ^ n) = n * (u ^ (n-1))* du_dx */
809 du_dx = Derivative(left,wrt,filter);
810 e1 = MakePower(left,CreateDTermInteger(ivalue-1));
811 e1 = MakeMultiply(CreateDTermInteger(ivalue),MakeMultiply(e1,du_dx));
812 return e1;
813 }
814 }
815 /* d/dx (u ^ n) = n * (u ^ (n-1))* du_dx */
816 if (right->t==e_real) {
817 rvalue = R_TERM(right)->value;
818 du_dx = Derivative(left,wrt,filter);
819 e1 = MakePower(left,CreateDTermReal(rvalue-1.0,NULL));
820 e1 = MakeMultiply(CreateDTermReal(rvalue,NULL),MakeMultiply(e1,du_dx));
821 return e1;
822 }
823
824 du_dx = Derivative(left,wrt,filter);
825 dv_dx = Derivative(right,wrt,filter);
826 e0 = MakeFunc(F_LN,TermBinRight(term));
827 e0 = MakeMultiply(dv_dx,e0);
828 e1 = MakeDivide(TermBinRight(term),TermBinRight(term));
829 e1 = MakeMultiply(du_dx,e1);
830 e0 = MakeAdd(e0,e1);
831 e0 = MakeMultiply(e0,MakePower(TermBinRight(term),TermBinRight(term)));
832 return e0;
833 }
834
835 case e_func: { /* deal with all the recognized functions */
836 enum Func_enum id;
837 Term *left;
838 Term *e;
839
840 left = CloneDTerm(TermFuncLeft(term));
841 du_dx = Derivative(left,wrt,filter);
842 switch(id=FuncId(TermFunc(term))){ /* d (sin(u))/dx = cos(u).du_dx */
843 case F_SIN:
844 e = MakeFunc(F_COS,left);
845 e = MakeMultiply(du_dx,e);
846 return e;
847 case F_SINH:
848 e = MakeFunc(F_SINH,left);
849 e = MakeMultiply(du_dx,e);
850 return e;
851 case F_COS:
852 e = MakeFunc(F_SIN,left);
853 e = MakeMultiply(du_dx,e);
854 e = MakeNegation(e);
855 return e;
856 case F_COSH:
857 e = MakeFunc(F_SINH,left);
858 e = MakeMultiply(du_dx,e);
859 return e;
860 case F_TAN:
861 e = MakeFunc(F_COS,left);
862 e = MakeFunc(F_SQR,left);
863 e = MakeDivide(du_dx,e);
864 return e;
865 case F_TANH:
866 e = MakeFunc(F_COSH,left);
867 e = MakeFunc(F_SQR,e);
868 e = MakeDivide(du_dx,e);
869 return e;
870 case F_ARCSIN:
871 e = MakeFunc(F_SQR,left);
872 e = MakeSubtract(K_TERM(&One),e);
873 e = MakeFunc(F_SQRT,e);
874 e = MakeDivide(du_dx,e);
875 return e;
876 case F_ARCSINH:
877 e = MakeFunc(F_SQR,left);
878 e = MakeAdd(e,K_TERM(&One));
879 e = MakeFunc(F_SQRT,e);
880 e = MakeDivide(du_dx,e);
881 return e;
882 case F_ARCCOS:
883 e = MakeFunc(F_SQR,left);
884 e = MakeSubtract(K_TERM(&One),e);
885 e = MakeFunc(F_SQRT,e);
886 e = MakeDivide(du_dx,e);
887 e = MakeNegation(e);
888 return e;
889 case F_ARCCOSH:
890 e = MakeFunc(F_SQR,left);
891 e = MakeSubtract(e,K_TERM(&One));
892 e = MakeFunc(F_SQRT,e);
893 e = MakeDivide(du_dx,e);
894 return e;
895 case F_ARCTAN:
896 e = MakeFunc(F_SQR,left);
897 e = MakeAdd(K_TERM(&One),e);
898 e = MakeDivide(du_dx,e);
899 return e;
900 case F_ARCTANH:
901 e = MakeFunc(F_SQR,left);
902 e = MakeSubtract(K_TERM(&One),e);
903 e = MakeDivide(du_dx,e);
904 return e;
905 #ifdef HAVE_ERF
906 case F_ERF:
907 e = MakeFunc(F_SQR,left);
908 e = MakeNegation(e);
909 e = MakeFunc(F_EXP,e);
910 e = MakeMultiply(e,CreateDTermReal(F_ERF_COEF,NULL));
911 e = MakeMultiply(du_dx,e);
912 return e;
913 #endif /* HAVE_ERF */
914 case F_EXP:
915 e = MakeMultiply(du_dx,term);
916 return e;
917 case F_LN:
918 e = MakeDivide(du_dx,left);
919 return e;
920 case F_LOG10:
921 e = MakeDivide(du_dx,left);
922 e = MakeMultiply(e,CreateDTermReal(F_LOG10_COEF,NULL));
923 return e;
924 case F_SQR:
925 e = MakeMultiply(K_TERM(&Two),left);
926 e = MakeMultiply(du_dx,e); /* for prettyness should swap */
927 return e;
928 case F_SQRT: /* d(u ^ (1/2))/dx = du_dx / ( 2 * u^(1/2)) */
929 e = MakeMultiply(K_TERM(&Two),term);
930 e = MakeDivide(du_dx,e); /* note that term is used */
931 return e;
932 case F_CUBE: /* d(u^3)/dx = 3 . u^2. du_dx */
933 e = MakeFunc(F_SQR,left);
934 e = MakeMultiply(K_TERM(&Three),e);
935 return e;
936 /*
937 * more efficient that using sqr(u)^1/3, as other
938 * simplifications become evident later.
939 */
940 case F_CBRT: /* d(u ^ (1/3))/dx = du_dx / (3 * u^(2/3)) */
941 e = CreateDTermReal(2.0/3.0,NULL);
942 e = MakePower(left,e);
943 e = MakeMultiply(K_TERM(&Three),e);
944 e = MakeDivide(du_dx,e);
945 return e;
946 default:
947 FPRINTF(ASCERR,"Unknown function type in case e_func\n");
948 return NULL;
949 }
950 } /* end of e_func case */
951
952 case e_uminus: { /* d (-u)/dx = - du_dx */
953 Term *left;
954 left = CloneDTerm(TermUniLeft(term));
955 du_dx = Derivative(left,wrt,filter);
956
957 if (du_dx==K_TERM(&Zero))
958 return K_TERM(&Zero);
959 return MakeNegation(du_dx);
960 }
961
962 default:
963 FPRINTF(ASCERR,"Unknown term type in Derivatives\n");
964 return NULL; /* only time that NULL should appear */
965 }
966 }
967
968 /*
969 * This function assumes that the derivative free store exists.
970 * It will only ReInit the store before starting. For the time
971 * In the final implementation this will be a pointer to an Instance,
972 * or a index into a table of instances. We allow the reset flag
973 * as we might choose not to reset after performing the derivatives on
974 * say the lhs of a relation, but wait until the rhs has been
975 * completed as well.
976 */
977
978 Term *TermDerivative(Term *term, unsigned long wrt,
979 int (*filter)(struct Instance *))
980 {
981 Term *result = NULL;
982 assert(term!=NULL);
983 result = Derivative(CloneDTerm(term),wrt,filter);
984 return result;
985 }
986
987 RelationINF *RelDerivative(RelationINF *rel,unsigned long wrt,
988 int (*filter)(struct Instance *))
989 {
990 Term *side; /* taken from the freestore */
991 RelationINF *result;
992 unsigned long len;
993
994 result = CreateRelInfix(NULL,rel->share->s.relop,NULL);
995
996 len = NumberVariables(rel);
997 if (!len) { /* no variables present */
998 RTOKEN(result).lhs_term = CopyDTerm(K_TERM(&Zero));
999 RTOKEN(result).rhs_term = CopyDTerm(K_TERM(&Zero));
1000 return result;
1001 }
1002 /*
1003 * 1) Initialize.
1004 * 2) From a = b , create a - b = 0; node allocated from freestore.
1005 * 3) Then differentiate this single function.
1006 * 4) Make a fresh copy for the user.
1007 * 5) Fix the varlist, and return the relation -- The user owns all.
1008 */
1009 FreeStore_ReInit(FreeStore_GetFreeStore());
1010 side = CreateDTermBinary(Infix_LhsSide(rel),
1011 e_minus,
1012 Infix_RhsSide(rel));
1013 RTOKEN(result).lhs_term = TermDerivative(side,wrt,filter);
1014 RTOKEN(result).lhs_term = CopyDTermTree(Infix_LhsSide(rel));
1015 RTOKEN(result).rhs_term = CopyDTerm(K_TERM(&Zero));
1016 result->vars = RedoVarList(result,rel->vars);
1017
1018 #ifdef FREESTORE_STATS
1019 FreeStore__Statistics(ASCERR,FreeStore_GetFreeStore());
1020 #endif
1021 return result;
1022 }
1023
1024 RelationINF *RelDeriveSloppy(RelationINF *rel,unsigned long wrt,
1025 int (*filter)(struct Instance *))
1026 {
1027 Term *side;
1028 RelationINF *result;
1029 unsigned long len;
1030
1031 result = CreateRelInfix(NULL,rel->share->s.relop,NULL);
1032 len = NumberVariables(rel);
1033 if (!len) { /* no variables present */
1034 RTOKEN(result).lhs_term = K_TERM(&Zero);
1035 RTOKEN(result).rhs_term = K_TERM(&Zero);
1036 return result;
1037 }
1038 /*
1039 * Initialize freestore
1040 * From a = b , create a - b = 0; node allocated from freestore
1041 * Then differentiate this single function.
1042 */
1043 FreeStore_ReInit(FreeStore_GetFreeStore());
1044 side = CreateDTermBinary(Infix_LhsSide(rel),
1045 e_minus,
1046 Infix_RhsSide(rel));
1047 RTOKEN(result).lhs_term = TermDerivative(side,wrt,filter);
1048 RTOKEN(result).rhs_term = K_TERM(&Zero);
1049 result->vars = gl_copy(rel->vars); /* not unique */
1050
1051 #ifdef FREESTORE_STATS
1052 FreeStore__Statistics(ASCERR,FreeStore_GetFreeStore());
1053 #endif
1054 return result;
1055 }
1056
1057 void RelDestroySloppy(RelationINF *rel)
1058 {
1059 Term *side;
1060 if (rel->vars) gl_destroy(rel->vars);
1061 side = A_TERM(FreeStoreCheckMem(UNION_TERM(Infix_LhsSide(rel))));
1062 side = A_TERM(FreeStoreCheckMem(UNION_TERM(Infix_RhsSide(rel))));
1063 }
1064
1065 void PrepareDerivatives(int setup,int n_buffers,int buffer_length)
1066 {
1067 static struct FreeStore *deriv_store = NULL;
1068 static struct FreeStore *previous_store = NULL;
1069
1070 if (setup) {
1071 previous_store = FreeStore_GetFreeStore();
1072 deriv_store = FreeStore_Create(n_buffers,buffer_length);
1073 FreeStore_SetFreeStore(deriv_store);
1074 InitConstantTerms();
1075 }
1076 else{ /* we are to shut down */
1077 if (deriv_store)
1078 FreeStore__BlastMem(deriv_store);
1079 FreeStore_SetFreeStore(previous_store);
1080 deriv_store = NULL;
1081 previous_store = NULL;
1082 }
1083 }
1084
1085
1086

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22