Annotation of OpenXM_contrib2/asir2000/builtin/pf.c, Revision 1.1.1.1
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/builtin/pf.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3: #include "math.h"
4: #include "parse.h"
5: #if 0
6: #include <alloca.h>
7: #endif
8:
9: double const_pi(),const_e();
10:
11: void make_ihyp(void);
12: void make_hyp(void);
13: void make_itri(void);
14: void make_tri(void);
15: void make_exp(void);
16: void simplify_pow(PFINS,Obj *);
17:
18: void Pfunctor(),Pargs(),Pfunargs(),Pvtype(),Pcall(),Pdeval();
19: void Pregister_handler();
20:
21: struct ftab puref_tab[] = {
22: {"functor",Pfunctor,1},
23: {"args",Pargs,1},
24: {"funargs",Pfunargs,1},
25: {"register_handler",Pregister_handler,1},
26: {"call",Pcall,2},
27: {"vtype",Pvtype,1},
28: {"deval",Pdeval,1},
29: {0,0,0},
30: };
31:
32: #if PARI
33: int p_pi(),p_e();
34: int p_log(),p_exp(),p_pow();
35: int p_sin(),p_cos(),p_tan(),p_asin(),p_acos(),p_atan();
36: int p_sinh(),p_cosh(),p_tanh(),p_asinh(),p_acosh(),p_atanh();
37: #else
38: int p_pi,p_e;
39: int p_log,p_exp,p_pow;
40: int p_sin,p_cos,p_tan,p_asin,p_acos,p_atan;
41: int p_sinh,p_cosh,p_tanh,p_asinh,p_acosh,p_atanh;
42: #endif
43:
44: static V *uarg,*darg;
45: static P x,y;
46: static PF pidef,edef;
47: static PF logdef,expdef,powdef;
48: static PF sindef,cosdef,tandef;
49: static PF asindef,acosdef,atandef;
50: static PF sinhdef,coshdef,tanhdef;
51: static PF asinhdef,acoshdef,atanhdef;
52:
53: #define OALLOC(p,n) ((p)=(Obj *)CALLOC((n),sizeof(Obj)))
54:
55: double const_pi() { return 3.14159265358979323846264338327950288; }
56: double const_e() { return 2.718281828459045235360287471352662497; }
57:
58: void pf_init() {
59: uarg = (V *)CALLOC(1,sizeof(V));
60: uarg[0] = &oVAR[26]; MKV(uarg[0],x);
61:
62: darg = (V *)CALLOC(2,sizeof(V));
63: darg[0] = &oVAR[26];
64: darg[1] = &oVAR[27]; MKV(darg[1],y);
65:
66: mkpf("@pi",0,0,0,(int (*)())p_pi,const_pi,0,&pidef);
67: mkpf("@e",0,0,0,(int (*)())p_e,const_e,0,&edef);
68:
69: mkpf("log",0,1,uarg,(int (*)())p_log,log,0,&logdef);
70: mkpf("exp",0,1,uarg,(int (*)())p_exp,exp,0,&expdef);
71: mkpf("pow",0,2,darg,(int (*)())p_pow,pow,(int (*)())simplify_pow,&powdef);
72:
73: mkpf("sin",0,1,uarg,(int (*)())p_sin,sin,0,&sindef);
74: mkpf("cos",0,1,uarg,(int (*)())p_cos,cos,0,&cosdef);
75: mkpf("tan",0,1,uarg,(int (*)())p_tan,tan,0,&tandef);
76: mkpf("asin",0,1,uarg,(int (*)())p_asin,asin,0,&asindef);
77: mkpf("acos",0,1,uarg,(int (*)())p_acos,acos,0,&acosdef);
78: mkpf("atan",0,1,uarg,(int (*)())p_atan,atan,0,&atandef);
79:
80: mkpf("sinh",0,1,uarg,(int (*)())p_sinh,sinh,0,&sinhdef);
81: mkpf("cosh",0,1,uarg,(int (*)())p_cosh,cosh,0,&coshdef);
82: mkpf("tanh",0,1,uarg,(int (*)())p_tanh,tanh,0,&tanhdef);
83: #if !defined(VISUAL)
84: mkpf("asinh",0,1,uarg,(int (*)())p_asinh,asinh,0,&asinhdef);
85: mkpf("acosh",0,1,uarg,(int (*)())p_acosh,acosh,0,&acoshdef);
86: mkpf("atanh",0,1,uarg,(int (*)())p_atanh,atanh,0,&atanhdef);
87: #endif
88: make_exp();
89: make_tri();
90: make_itri();
91: make_hyp();
92: #if !defined(VISUAL)
93: make_ihyp();
94: #endif
95: }
96:
97: void make_exp() {
98: V v;
99: P u,vexp,vlog,vpow;
100: Obj *args;
101:
102: mkpfins(expdef,uarg,&v); MKV(v,vexp);
103: mkpfins(powdef,darg,&v); MKV(v,vpow);
104: mkpfins(logdef,uarg,&v); MKV(v,vlog);
105:
106: /* d/dx(log(x)) = 1/x */
107: OALLOC(logdef->deriv,1); divr(CO,(Obj)ONE,(Obj)x,&logdef->deriv[0]);
108:
109: /* d/dx(exp(x)) = exp(x) */
110: OALLOC(expdef->deriv,1); expdef->deriv[0] = (Obj)vexp;
111:
112: /* d/dy(x^y) = log(x)*x^y */
113: OALLOC(powdef->deriv,2); mulp(CO,vpow,vlog,(P *)&powdef->deriv[1]);
114:
115: /* d/dx(x^y) = y*x^(y-1) */
116: args = (Obj *)ALLOCA(2*sizeof(Obj));
117: args[0] = (Obj)x; subp(CO,y,(P)ONE,(P *)&args[1]);
118: _mkpfins(powdef,args,&v); MKV(v,u);
119: mulr(CO,(Obj)u,(Obj)y,&powdef->deriv[0]);
120: }
121:
122: void make_tri() {
123: V v;
124: P vcos,vsin,vtan,t;
125:
126: mkpfins(cosdef,uarg,&v); MKV(v,vcos);
127: mkpfins(sindef,uarg,&v); MKV(v,vsin);
128: mkpfins(tandef,uarg,&v); MKV(v,vtan);
129:
130: /* d/dx(sin(x)) = cos(x) */
131: OALLOC(sindef->deriv,1); sindef->deriv[0] = (Obj)vcos;
132:
133: /* d/dx(cos(x)) = -sin(x) */
134: OALLOC(cosdef->deriv,1); chsgnp(vsin,(P *)&cosdef->deriv[0]);
135:
136: /* d/dx(tan(x)) = 1+tan(x)^2 */
137: OALLOC(tandef->deriv,1);
138: mulr(CO,(Obj)vtan,(Obj)vtan,(Obj *)&t); addp(CO,(P)ONE,t,(P *)&tandef->deriv[0]);
139: }
140:
141: void make_itri() {
142: P t,xx;
143: Q mtwo;
144: V v;
145: Obj *args;
146:
147: /* d/dx(asin(x)) = (1-x^2)^(-1/2) */
148: OALLOC(asindef->deriv,1);
149: args = (Obj *)ALLOCA(2*sizeof(Obj));
150: mulp(CO,x,x,&xx); subp(CO,(P)ONE,xx,(P *)&args[0]);
151: STOQ(-2,mtwo); divq(ONE,mtwo,(Q *)&args[1]);
152: _mkpfins(powdef,args,&v); MKV(v,t);
153: asindef->deriv[0] = (Obj)t;
154:
155: /* d/dx(acos(x)) = -(1-x^2)^(-1/2) */
156: OALLOC(acosdef->deriv,1); chsgnp((P)asindef->deriv[0],(P *)&acosdef->deriv[0]);
157:
158: /* d/dx(atan(x)) = 1/(x^2+1) */
159: OALLOC(atandef->deriv,1);
160: addp(CO,(P)ONE,xx,&t); divr(CO,(Obj)ONE,(Obj)t,&atandef->deriv[0]);
161: }
162:
163: void make_hyp() {
164: V v;
165: P vcosh,vsinh,vtanh,t;
166:
167: mkpfins(coshdef,uarg,&v); MKV(v,vcosh);
168: mkpfins(sinhdef,uarg,&v); MKV(v,vsinh);
169: mkpfins(tanhdef,uarg,&v); MKV(v,vtanh);
170:
171: /* d/dx(sinh(x)) = cosh(x) */
172: OALLOC(sinhdef->deriv,1); sinhdef->deriv[0] = (Obj)vcosh;
173:
174: /* d/dx(cosh(x)) = sinh(x) */
175: OALLOC(coshdef->deriv,1); coshdef->deriv[0] = (Obj)vsinh;
176:
177: /* d/dx(tanh(x)) = 1-tanh(x)^2 */
178: OALLOC(tanhdef->deriv,1);
179: mulr(CO,(Obj)vtanh,(Obj)vtanh,(Obj *)&t); subp(CO,(P)ONE,t,(P *)&tanhdef->deriv[0]);
180: }
181:
182: void make_ihyp() {
183: P t,xx;
184: Q mtwo;
185: V v;
186: Obj *args;
187:
188: /* d/dx(asinh(x)) = (1+x^2)^(-1/2) */
189: OALLOC(asinhdef->deriv,1);
190: args = (Obj *)ALLOCA(2*sizeof(Obj));
191: mulp(CO,x,x,&xx); addp(CO,(P)ONE,xx,(P *)&args[0]);
192: STOQ(-2,mtwo); divq(ONE,mtwo,(Q *)&args[1]);
193: _mkpfins(powdef,args,&v); MKV(v,t);
194: asinhdef->deriv[0] = (Obj)t;
195:
196: /* d/dx(acosh(x)) = (x^2-1)^(-1/2) */
197: OALLOC(acoshdef->deriv,1);
198: subp(CO,xx,(P)ONE,(P *)&args[0]);
199: _mkpfins(powdef,args,&v); MKV(v,t);
200: acoshdef->deriv[0] = (Obj)t;
201:
202: /* d/dx(atanh(x)) = 1/(1-x^2) */
203: OALLOC(atanhdef->deriv,1);
204: subp(CO,(P)ONE,xx,&t); divr(CO,(Obj)ONE,(Obj)t,&atanhdef->deriv[0]);
205: }
206:
207: void mkpow(vl,a,e,r)
208: VL vl;
209: Obj a;
210: Obj e;
211: Obj *r;
212: {
213: PFINS ins;
214: PFAD ad;
215:
216: ins = (PFINS)CALLOC(1,sizeof(PF)+2*sizeof(struct oPFAD));
217: ins->pf = powdef; ad = ins->ad;
218: ad[0].d = 0; ad[0].arg = a; ad[1].d = 0; ad[1].arg = e;
219: simplify_ins(ins,r);
220: }
221:
222: void simplify_pow(ins,rp)
223: PFINS ins;
224: Obj *rp;
225: {
226: PF pf;
227: PFAD ad;
228: Obj a0,a1;
229: V v;
230: P t;
231:
232: pf = ins->pf; ad = ins->ad; a0 = ad[0].arg; a1 = ad[1].arg;
233: if ( !a1 )
234: *rp = (Obj)ONE;
235: else if ( !a0 )
236: *rp = 0;
237: else if ( NUM(a1) && INT(a1) )
238: arf_pwr(CO,a0,a1,rp);
239: else {
240: instov(ins,&v); MKV(v,t); *rp = (Obj)t;
241: }
242: }
243:
244: #define ISPFINS(p)\
245: (p)&&(ID(p) == O_P)&&((int)VR((P)p)->attr!=V_PF)&&\
246: UNIQ(DEG(DC((P)p)))&&UNIQ(COEF(DC((P)p)))
247:
248: void Pfunctor(arg,rp)
249: NODE arg;
250: P *rp;
251: {
252: P p;
253: FUNC t;
254: PF pf;
255: PFINS ins;
256:
257: p = (P)ARG0(arg);
258: if ( !ISPFINS(p) )
259: *rp = 0;
260: else {
261: ins = (PFINS)VR(p)->priv; pf = ins->pf;
262: t = (FUNC)MALLOC(sizeof(struct oFUNC));
263: t->name = pf->name; t->id = A_PURE; t->argc = pf->argc;
264: t->f.puref = pf;
265: makesrvar(t,rp);
266: }
267: }
268:
269: void Pargs(arg,rp)
270: NODE arg;
271: LIST *rp;
272: {
273: P p;
274: PF pf;
275: PFAD ad;
276: PFINS ins;
277: NODE n,n0;
278: int i;
279:
280: p = (P)ARG0(arg);
281: if ( !ISPFINS(p) )
282: *rp = 0;
283: else {
284: ins = (PFINS)VR(p)->priv; ad = ins->ad; pf = ins->pf;
285: for ( i = 0, n0 = 0; i < pf->argc; i++ ) {
286: NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
287: }
288: if ( n0 )
289: NEXT(n) = 0;
290: MKLIST(*rp,n0);
291: }
292: }
293:
294: void Pfunargs(arg,rp)
295: NODE arg;
296: LIST *rp;
297: {
298: P p;
299: P f;
300: FUNC t;
301: PF pf;
302: PFINS ins;
303: PFAD ad;
304: NODE n,n0;
305: int i;
306:
307: p = (P)ARG0(arg);
308: if ( !ISPFINS(p) )
309: *rp = 0;
310: else {
311: ins = (PFINS)VR(p)->priv; ad = ins->ad; pf = ins->pf;
312: t = (FUNC)MALLOC(sizeof(struct oFUNC));
313: t->name = pf->name; t->id = A_PURE; t->argc = pf->argc;
314: t->f.puref = pf;
315: makesrvar(t,&f);
316: n0 = 0; NEXTNODE(n0,n); BDY(n) = (pointer)f;
317: for ( i = 0; i < pf->argc; i++ ) {
318: NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
319: }
320: NEXT(n) = 0;
321: MKLIST(*rp,n0);
322: }
323: }
324:
325: void Pvtype(arg,rp)
326: NODE arg;
327: Q *rp;
328: {
329: P p;
330:
331: p = (P)ARG0(arg);
332: if ( !p || ID(p) != O_P )
333: *rp = 0;
334: else
335: STOQ((int)VR(p)->attr,*rp);
336: }
337:
338: extern FUNC registered_handler;
339:
340: void Pregister_handler(arg,rp)
341: NODE arg;
342: Q *rp;
343: {
344: P p;
345: V v;
346: FUNC func;
347:
348: p = (P)ARG0(arg);
349: if ( !p )
350: registered_handler = 0;
351: else if ( OID(p) != 2 )
352: error("register_hanlder : invalid argument");
353: v = VR(p);
354: if ( (int)v->attr != V_SR )
355: error("register_hanlder : no such function");
356: else {
357: func = (FUNC)v->priv;
358: if ( func->argc )
359: error("register_hanlder : the function must be with no argument");
360: else {
361: registered_handler = func;
362: *rp = ONE;
363: }
364: }
365: }
366:
367: void Pcall(arg,rp)
368: NODE arg;
369: Obj *rp;
370: {
371: P p;
372: V v;
373:
374: p = (P)ARG0(arg);
375: if ( !p || OID(p) != 2 )
376: error("call : invalid argument");
377: v = VR(p);
378: if ( (int)v->attr != V_SR )
379: error("call : no such function");
380:
381: else
382: *rp = (Obj)bevalf((FUNC)v->priv,BDY((LIST)ARG1(arg)));
383: }
384:
385: void Pdeval(arg,rp)
386: NODE arg;
387: Obj *rp;
388: {
389: asir_assert(ARG0(arg),O_R,"deval");
390: devalr(CO,(Obj)ARG0(arg),rp);
391: }
392:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>