![[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.6 ! ohara 1: /* $OpenXM: OpenXM/src/ox_gsl/ox_eval.c,v 1.5 2018/04/13 16:51:42 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) {
121: e->name = s;
122: e->v = v;
123: e->f = NULL;
124: local_dic_counter++;
125: return 1;
126: }
127: return 0;
128: }
129:
130: void init_dic()
131: {
132: int i;
133: for(i=0; i<local_dic_counter; i++) {
134: free((local_dic+i)->name);
135: }
136: local_dic_counter=0;
137: memset(local_dic, 0, sizeof(entry)*LOCAL_DIC_SIZE);
138: }
139:
1.3 ohara 140: void replace(int n, ...)
141: {
142: char *s;
143: double d;
144: va_list ap;
145: va_start(ap,n);
146: for(init_dic(); n>0; n--) {
147: s = va_arg(ap, char *);
148: d = va_arg(ap, double);
149: register_entry(s,d);
150: }
151: va_end(ap);
152: }
153:
154: void replace2(int n, char *s[], double v[])
155: {
156: int i;
157: init_dic();
158: for(i=0; i<n; i++) {
159: register_entry(s[i],v[i]);
160: }
161: }
162:
1.4 ohara 163: static entry *find_entry(cmo *node, int len, entry *dic)
1.1 ohara 164: {
165: char *s;
166: entry *e;
167: if(node->tag == CMO_STRING) {
168: s = ((cmo_string *)node)->s;
169: }else if(node->tag == CMO_INDETERMINATE) {
170: s = cmo_indeterminate_get_name((cmo_indeterminate *)node);
171: }else {
172: return NULL;
173: }
174: for(e=dic; e->name != NULL; e++) {
1.4 ohara 175: if(strcmp(e->name,s)==0 && (len<0 || len==e->n_args)) {
1.1 ohara 176: return e;
177: }
178: }
179: return NULL;
180: }
181:
182: static int entry_function(entry *e, cmo_list *args, double *retval)
183: {
184: int i,m,n;
185: double *dargs;
186: n=e->n_args;
187: if(n<0 || n>list_length(args)) {
188: /* arguments are mismatched */
189: return 0;
190: }
191: if(n>5) {
192: /* too many arguments */
193: return 0;
194: }
195: dargs = (double *)alloca(n*sizeof(double));
196: for(i=0; i<n; i++) {
197: if(eval_cmo(list_nth(args, i), &dargs[i])==0) {
198: return 0;
199: }
200: }
201: switch(n) {
202: case 0:
203: *retval = e->f();
204: break;
205: case 1:
206: *retval = e->f(dargs[0]);
207: break;
208: case 2:
209: *retval = e->f(dargs[0],dargs[1]);
210: break;
211: case 3:
212: *retval = e->f(dargs[0],dargs[1],dargs[2]);
213: break;
214: case 4:
215: *retval = e->f(dargs[0],dargs[1],dargs[2],dargs[3]);
216: break;
217: case 5:
218: *retval = e->f(dargs[0],dargs[1],dargs[2],dargs[3],dargs[4]);
219: break;
220: default:
221: return 0;
222: }
223: return 1;
224: }
225:
226: int entry_value(entry *e, double *retval)
227: {
228: *retval = e->v;
229: return 1;
230: }
231:
232: /*
233: 返り値: 評価が成功したか、 1: 成功, 0: 失敗
234: 評価された値: *retval に格納
235: */
236: static int eval_cmo_tree(cmo_tree* t, double *retval)
237: {
1.4 ohara 238: entry *e = find_entry((cmo *)t->name,list_length(t->leaves),global_dic);
239: if (FUNCTION_P(e)) {
240: return entry_function(e,t->leaves,retval);
241: }else if (VALUE_P(e)) {
242: return entry_value(e, retval);
243: }
1.1 ohara 244: return 0;
245: }
246:
247: static int eval_cmo_indeterminate(cmo_indeterminate *c, double *retval)
248: {
1.4 ohara 249: entry *e = find_entry((cmo *)c,-1,local_dic);
1.1 ohara 250: if (VALUE_P(e)) {
251: return entry_value(e,retval);
252: }
253: return 0;
254: }
255:
1.2 ohara 256: static double mypow(double x, int n)
257: {
258: int i,k,f=0;
259: double s,t;
260: if (n==0) {
261: return 1;
262: }else if (n<0) {
263: n=-n;
264: f=1;
265: }
266: /* t=x^(2^i) */
267: s=1;
268: t=x;
269: for(k=n; k!=0; k=k>>1) {
270: if(k&1) {
271: s*=t;
272: }
273: t=t*t;
274: }
275: if (f>0) {
276: s = 1/s;
277: }
278: return s;
279: }
280:
281: static int eval_cmo_polynomial_in_one_variable(cmo_polynomial_in_one_variable* c, double vars[], int n, double *retval)
282: {
283: int i;
284: cell *cc;
285: double r,s=0;
286: double x = vars[c->var];
287: double *d=(double *)calloc(c->length, sizeof(double));
288: for(i=0; i<c->length; i++) {
289: cc = list_nth_cell((cmo_list *)c, i);
290: if (cc->cmo->tag == CMO_POLYNOMIAL_IN_ONE_VARIABLE) {
291: if(eval_cmo_polynomial_in_one_variable((cmo_polynomial_in_one_variable*)cc->cmo,vars,n,&r)==0) {
292: return 0;
293: }
294: }else {
295: if(eval_cmo(cc->cmo,&r)==0) {
296: return 0;
297: }
298: }
299: s += mypow(x,cc->exp)*r;
300: }
301: *retval = s;
302: return 1;
303: }
304:
305: static int eval_cmo_recursive_polynomial(cmo_recursive_polynomial* c, double *retval)
306: {
307: int i,n;
308: double *vars;
309: entry *e;
310: switch(c->coef->tag) {
311: case CMO_POLYNOMIAL_IN_ONE_VARIABLE:
312: n=list_length(c->ringdef);
313: if(local_dic_counter<n) {
314: return 0; /* 自由変数が残る */
315: }
316: vars=(double *)calloc(n,sizeof(double));
317: for(i=0; i<n; i++) {
1.4 ohara 318: e = find_entry(list_nth(c->ringdef,i),-1,local_dic);
1.2 ohara 319: if(e == NULL) {
320: free(vars);
321: return 0; /* failed */
322: }
323: entry_value(e, &vars[i]);
324: }
325: return eval_cmo_polynomial_in_one_variable((cmo_polynomial_in_one_variable*)c->coef,vars,n,retval);
326: case CMO_DISTRIBUTED_POLYNOMIAL:
327: return 0; /* failed */
328: default: /* cmo_zz, cmo_qq, cmo_double, ... */
329: return eval_cmo(c->coef,retval);
330: }
331: }
332:
1.1 ohara 333: int eval_cmo(cmo *c, double *retval)
334: {
335: int tag = c->tag;
336: switch(c->tag) {
1.5 ohara 337: case CMO_NULL:
1.1 ohara 338: case CMO_ZERO:
339: *retval = 0;
340: break;
341: case CMO_INT32:
342: *retval = (double)((cmo_int32 *)c)->i;
343: break;
344: case CMO_ZZ:
345: *retval = mpz_get_d(((cmo_zz *)c)->mpz);
346: break;
347: case CMO_QQ:
348: *retval = mpq_get_d(((cmo_qq *)c)->mpq);
349: break;
350: case CMO_IEEE_DOUBLE_FLOAT:
351: case CMO_64BIT_MACHINE_DOUBLE:
352: *retval = ((cmo_double *)c)->d;
353: break;
354: case CMO_BIGFLOAT:
355: *retval = mpfr_get_d(((cmo_bf *)c)->mpfr,GMP_RNDN);
356: break;
357: case CMO_TREE:
358: return eval_cmo_tree((cmo_tree *)c,retval);
359: break;
360: case CMO_INDETERMINATE:
361: return eval_cmo_indeterminate((cmo_indeterminate *)c,retval);
362: break;
1.2 ohara 363: case CMO_RECURSIVE_POLYNOMIAL:
364: return eval_cmo_recursive_polynomial((cmo_recursive_polynomial *)c,retval);
365: break;
1.1 ohara 366: default:
367: /* 変換できない型 */
368: return 0; /* error */
369: }
370: return 1;
371: }