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