![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ox_eval.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / ox_gsl|
Return to ox_eval.c CVS log | Up to [local] / OpenXM / src / ox_gsl |
1.7 ! ohara 1: /* $OpenXM: OpenXM/src/ox_gsl/ox_eval.c,v 1.6 2018/04/17 02:50:07 ohara Exp $ */
1.1 ohara 2:
3: #include <stdio.h>
4: #include <stdlib.h>
1.3 ohara 5: #include <stdarg.h>
1.1 ohara 6: #include <string.h>
7: #include <math.h>
8: #include "ox_toolkit.h"
9:
10: /*
11: Usage:
12:
13: double d;
1.3 ohara 14: replace(3,"x",1.25,"y",-2.0, "z", 2.1);
1.1 ohara 15: if(eval_cmo(your_cmo_tree,&d)==0) goto_error();
16: */
17:
18: #define FUNCTION_P(e) (((e)!=NULL) && ((e)->f != NULL))
19: #define VALUE_P(e) (((e)!=NULL) && ((e)->f == NULL))
20:
1.2 ohara 21: #define FAILED 0
22: #define SUCCEED 1
23:
1.3 ohara 24: void replace(int n, ...);
25: void replace2(int n, char *s[], double v[]);
1.1 ohara 26: int eval_cmo(cmo *c, double *retval);
27:
28: static double op_add(double x,double y)
29: {
30: return x+y;
31: }
32:
33: static double op_sub(double x,double y)
34: {
35: return x-y;
36: }
37:
38: static double op_mul(double x,double y)
39: {
40: return x*y;
41: }
42:
43: static double op_div(double x,double y)
44: {
45: return x/y;
46: }
47:
48: static double op_negative(double x)
49: {
50: return -x;
51: }
52:
1.6 ohara 53: static double op_parentheses(double x)
54: {
55: return x;
56: }
57:
1.1 ohara 58: /* 定数は引数なしの関数として実現する。*/
59: typedef struct {
60: char *name;
61: double v;
62: double (*f)();
63: int n_args; /* number of args */
64: } entry;
65:
66: /*
67: グローバル辞書: sin(x), cos(x), exp(x), lgamma(x) などの関数や、pi,e といった定数を扱うための辞書
68: (グローバル変数で保持、システム固有)
69:
70: ローカル辞書: [x -> 1/2, y-> 0.5] などの置き換えのため。
71: */
72: entry global_dic[512] = {
73: {"sin",0,sin,1},
74: {"cos",0,cos,1},
1.4 ohara 75: {"tan",0,tan,1},
76: {"sinh",0,sinh,1},
77: {"cosh",0,cosh,1},
78: {"tanh",0,tanh,1},
79: {"asin",0,asin,1},
80: {"acos",0,acos,1},
81: {"asinh",0,asinh,1},
82: {"acosh",0,acosh,1},
83: {"erf",0,erf,1},
1.1 ohara 84: {"exp",0,exp,1},
1.4 ohara 85: {"exp2",0,exp2,1},
1.1 ohara 86: {"log",0,log,1},
1.4 ohara 87: {"log2",0,log2,1},
88: {"log10",0,log10,1},
89: {"gamma",0,gamma,1},
90: {"lgamma",0,lgamma,1},
91: {"sqrt",0,sqrt,1},
92: {"cbrt",0,cbrt,1},
93: {"fabs",0,fabs,1},
94: {"j0",0,j0,1},
95: {"j1",0,j1,1},
96: {"y0",0,y0,1},
97: {"y1",0,y1,1},
1.6 ohara 98: {"()", 0,op_parentheses,1},
1.4 ohara 99: {"-", 0,op_negative,1},
1.1 ohara 100: {"+", 0,op_add,2},
101: {"-", 0,op_sub,2},
102: {"*", 0,op_mul,2},
103: {"/", 0,op_div,2},
104: {"^", 0,pow,2},
1.4 ohara 105: {"pow", 0,pow,2},
106: {"@e", M_E, NULL,0},
107: {"@pi", M_PI,NULL,0},
1.1 ohara 108: {NULL,0,NULL,0}
109: };
110:
111: #define LOCAL_DIC_SIZE 256
112: static entry local_dic[LOCAL_DIC_SIZE] = {
113: {NULL,0,NULL,0}
114: };
115: static int local_dic_counter;
116:
117: int register_entry(char *s, double v)
118: {
119: entry *e = &local_dic[local_dic_counter];
120: if(local_dic_counter<LOCAL_DIC_SIZE-1) {
1.7 ! ohara 121: e->name = (char *)malloc(strlen(s)+1);
! 122: strcpy(e->name, s);
1.1 ohara 123: e->v = v;
1.7 ! ohara 124: e->f = NULL;
1.1 ohara 125: local_dic_counter++;
126: return 1;
127: }
128: return 0;
129: }
130:
131: void init_dic()
132: {
133: int i;
134: for(i=0; i<local_dic_counter; i++) {
135: free((local_dic+i)->name);
136: }
137: local_dic_counter=0;
138: memset(local_dic, 0, sizeof(entry)*LOCAL_DIC_SIZE);
139: }
140:
1.3 ohara 141: void replace(int n, ...)
142: {
143: char *s;
144: double d;
145: va_list ap;
146: va_start(ap,n);
147: for(init_dic(); n>0; n--) {
148: s = va_arg(ap, char *);
149: d = va_arg(ap, double);
150: register_entry(s,d);
151: }
152: va_end(ap);
153: }
154:
155: void replace2(int n, char *s[], double v[])
156: {
157: int i;
158: init_dic();
159: for(i=0; i<n; i++) {
160: register_entry(s[i],v[i]);
161: }
162: }
163:
1.4 ohara 164: static entry *find_entry(cmo *node, int len, entry *dic)
1.1 ohara 165: {
166: char *s;
167: entry *e;
168: if(node->tag == CMO_STRING) {
169: s = ((cmo_string *)node)->s;
170: }else if(node->tag == CMO_INDETERMINATE) {
171: s = cmo_indeterminate_get_name((cmo_indeterminate *)node);
172: }else {
173: return NULL;
174: }
175: for(e=dic; e->name != NULL; e++) {
1.4 ohara 176: if(strcmp(e->name,s)==0 && (len<0 || len==e->n_args)) {
1.1 ohara 177: return e;
178: }
179: }
180: return NULL;
181: }
182:
183: static int entry_function(entry *e, cmo_list *args, double *retval)
184: {
185: int i,m,n;
186: double *dargs;
187: n=e->n_args;
188: if(n<0 || n>list_length(args)) {
189: /* arguments are mismatched */
190: return 0;
191: }
192: if(n>5) {
193: /* too many arguments */
194: return 0;
195: }
196: dargs = (double *)alloca(n*sizeof(double));
197: for(i=0; i<n; i++) {
198: if(eval_cmo(list_nth(args, i), &dargs[i])==0) {
199: return 0;
200: }
201: }
202: switch(n) {
203: case 0:
204: *retval = e->f();
205: break;
206: case 1:
207: *retval = e->f(dargs[0]);
208: break;
209: case 2:
210: *retval = e->f(dargs[0],dargs[1]);
211: break;
212: case 3:
213: *retval = e->f(dargs[0],dargs[1],dargs[2]);
214: break;
215: case 4:
216: *retval = e->f(dargs[0],dargs[1],dargs[2],dargs[3]);
217: break;
218: case 5:
219: *retval = e->f(dargs[0],dargs[1],dargs[2],dargs[3],dargs[4]);
220: break;
221: default:
222: return 0;
223: }
224: return 1;
225: }
226:
227: int entry_value(entry *e, double *retval)
228: {
229: *retval = e->v;
230: return 1;
231: }
232:
233: /*
234: 返り値: 評価が成功したか、 1: 成功, 0: 失敗
235: 評価された値: *retval に格納
236: */
237: static int eval_cmo_tree(cmo_tree* t, double *retval)
238: {
1.4 ohara 239: entry *e = find_entry((cmo *)t->name,list_length(t->leaves),global_dic);
240: if (FUNCTION_P(e)) {
241: return entry_function(e,t->leaves,retval);
242: }else if (VALUE_P(e)) {
243: return entry_value(e, retval);
244: }
1.1 ohara 245: return 0;
246: }
247:
248: static int eval_cmo_indeterminate(cmo_indeterminate *c, double *retval)
249: {
1.4 ohara 250: entry *e = find_entry((cmo *)c,-1,local_dic);
1.1 ohara 251: if (VALUE_P(e)) {
252: return entry_value(e,retval);
253: }
254: return 0;
255: }
256:
1.2 ohara 257: static double mypow(double x, int n)
258: {
259: int i,k,f=0;
260: double s,t;
261: if (n==0) {
262: return 1;
263: }else if (n<0) {
264: n=-n;
265: f=1;
266: }
267: /* t=x^(2^i) */
268: s=1;
269: t=x;
270: for(k=n; k!=0; k=k>>1) {
271: if(k&1) {
272: s*=t;
273: }
274: t=t*t;
275: }
276: if (f>0) {
277: s = 1/s;
278: }
279: return s;
280: }
281:
282: static int eval_cmo_polynomial_in_one_variable(cmo_polynomial_in_one_variable* c, double vars[], int n, double *retval)
283: {
284: int i;
285: cell *cc;
286: double r,s=0;
287: double x = vars[c->var];
288: double *d=(double *)calloc(c->length, sizeof(double));
289: for(i=0; i<c->length; i++) {
290: cc = list_nth_cell((cmo_list *)c, i);
291: if (cc->cmo->tag == CMO_POLYNOMIAL_IN_ONE_VARIABLE) {
292: if(eval_cmo_polynomial_in_one_variable((cmo_polynomial_in_one_variable*)cc->cmo,vars,n,&r)==0) {
293: return 0;
294: }
295: }else {
296: if(eval_cmo(cc->cmo,&r)==0) {
297: return 0;
298: }
299: }
300: s += mypow(x,cc->exp)*r;
301: }
302: *retval = s;
303: return 1;
304: }
305:
306: static int eval_cmo_recursive_polynomial(cmo_recursive_polynomial* c, double *retval)
307: {
308: int i,n;
309: double *vars;
310: entry *e;
311: switch(c->coef->tag) {
312: case CMO_POLYNOMIAL_IN_ONE_VARIABLE:
313: n=list_length(c->ringdef);
314: if(local_dic_counter<n) {
315: return 0; /* 自由変数が残る */
316: }
317: vars=(double *)calloc(n,sizeof(double));
318: for(i=0; i<n; i++) {
1.4 ohara 319: e = find_entry(list_nth(c->ringdef,i),-1,local_dic);
1.2 ohara 320: if(e == NULL) {
321: free(vars);
322: return 0; /* failed */
323: }
324: entry_value(e, &vars[i]);
325: }
326: return eval_cmo_polynomial_in_one_variable((cmo_polynomial_in_one_variable*)c->coef,vars,n,retval);
327: case CMO_DISTRIBUTED_POLYNOMIAL:
328: return 0; /* failed */
329: default: /* cmo_zz, cmo_qq, cmo_double, ... */
330: return eval_cmo(c->coef,retval);
331: }
332: }
333:
1.1 ohara 334: int eval_cmo(cmo *c, double *retval)
335: {
336: int tag = c->tag;
337: switch(c->tag) {
1.5 ohara 338: case CMO_NULL:
1.1 ohara 339: case CMO_ZERO:
340: *retval = 0;
341: break;
342: case CMO_INT32:
343: *retval = (double)((cmo_int32 *)c)->i;
344: break;
345: case CMO_ZZ:
346: *retval = mpz_get_d(((cmo_zz *)c)->mpz);
347: break;
348: case CMO_QQ:
349: *retval = mpq_get_d(((cmo_qq *)c)->mpq);
350: break;
351: case CMO_IEEE_DOUBLE_FLOAT:
352: case CMO_64BIT_MACHINE_DOUBLE:
353: *retval = ((cmo_double *)c)->d;
354: break;
355: case CMO_BIGFLOAT:
356: *retval = mpfr_get_d(((cmo_bf *)c)->mpfr,GMP_RNDN);
357: break;
358: case CMO_TREE:
359: return eval_cmo_tree((cmo_tree *)c,retval);
360: break;
361: case CMO_INDETERMINATE:
362: return eval_cmo_indeterminate((cmo_indeterminate *)c,retval);
363: break;
1.2 ohara 364: case CMO_RECURSIVE_POLYNOMIAL:
365: return eval_cmo_recursive_polynomial((cmo_recursive_polynomial *)c,retval);
366: break;
1.1 ohara 367: default:
368: /* 変換できない型 */
369: return 0; /* error */
370: }
371: return 1;
372: }