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