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