Annotation of OpenXM_contrib2/asir2000/builtin/pf.c, Revision 1.2
1.2 ! noro 1: /*
! 2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
! 3: * All rights reserved.
! 4: *
! 5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
! 6: * non-exclusive and royalty-free license to use, copy, modify and
! 7: * redistribute, solely for non-commercial and non-profit purposes, the
! 8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
! 9: * conditions of this Agreement. For the avoidance of doubt, you acquire
! 10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
! 11: * third party developer retains all rights, including but not limited to
! 12: * copyrights, in and to the SOFTWARE.
! 13: *
! 14: * (1) FLL does not grant you a license in any way for commercial
! 15: * purposes. You may use the SOFTWARE only for non-commercial and
! 16: * non-profit purposes only, such as academic, research and internal
! 17: * business use.
! 18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
! 19: * international copyright treaties. If you make copies of the SOFTWARE,
! 20: * with or without modification, as permitted hereunder, you shall affix
! 21: * to all such copies of the SOFTWARE the above copyright notice.
! 22: * (3) An explicit reference to this SOFTWARE and its copyright owner
! 23: * shall be made on your publication or presentation in any form of the
! 24: * results obtained by use of the SOFTWARE.
! 25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
! 26: * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
! 27: * for such modification or the source code of the modified part of the
! 28: * SOFTWARE.
! 29: *
! 30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
! 31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
! 32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
! 33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
! 34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
! 35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
! 36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
! 37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
! 38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
! 39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
! 40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
! 41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
! 42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
! 43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
! 44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
! 45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
! 46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
! 47: *
! 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/pf.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $
! 49: */
1.1 noro 50: #include "ca.h"
51: #include "math.h"
52: #include "parse.h"
53: #if 0
54: #include <alloca.h>
55: #endif
56:
57: double const_pi(),const_e();
58:
59: void make_ihyp(void);
60: void make_hyp(void);
61: void make_itri(void);
62: void make_tri(void);
63: void make_exp(void);
64: void simplify_pow(PFINS,Obj *);
65:
66: void Pfunctor(),Pargs(),Pfunargs(),Pvtype(),Pcall(),Pdeval();
67: void Pregister_handler();
68:
69: struct ftab puref_tab[] = {
70: {"functor",Pfunctor,1},
71: {"args",Pargs,1},
72: {"funargs",Pfunargs,1},
73: {"register_handler",Pregister_handler,1},
74: {"call",Pcall,2},
75: {"vtype",Pvtype,1},
76: {"deval",Pdeval,1},
77: {0,0,0},
78: };
79:
80: #if PARI
81: int p_pi(),p_e();
82: int p_log(),p_exp(),p_pow();
83: int p_sin(),p_cos(),p_tan(),p_asin(),p_acos(),p_atan();
84: int p_sinh(),p_cosh(),p_tanh(),p_asinh(),p_acosh(),p_atanh();
85: #else
86: int p_pi,p_e;
87: int p_log,p_exp,p_pow;
88: int p_sin,p_cos,p_tan,p_asin,p_acos,p_atan;
89: int p_sinh,p_cosh,p_tanh,p_asinh,p_acosh,p_atanh;
90: #endif
91:
92: static V *uarg,*darg;
93: static P x,y;
94: static PF pidef,edef;
95: static PF logdef,expdef,powdef;
96: static PF sindef,cosdef,tandef;
97: static PF asindef,acosdef,atandef;
98: static PF sinhdef,coshdef,tanhdef;
99: static PF asinhdef,acoshdef,atanhdef;
100:
101: #define OALLOC(p,n) ((p)=(Obj *)CALLOC((n),sizeof(Obj)))
102:
103: double const_pi() { return 3.14159265358979323846264338327950288; }
104: double const_e() { return 2.718281828459045235360287471352662497; }
105:
106: void pf_init() {
107: uarg = (V *)CALLOC(1,sizeof(V));
108: uarg[0] = &oVAR[26]; MKV(uarg[0],x);
109:
110: darg = (V *)CALLOC(2,sizeof(V));
111: darg[0] = &oVAR[26];
112: darg[1] = &oVAR[27]; MKV(darg[1],y);
113:
114: mkpf("@pi",0,0,0,(int (*)())p_pi,const_pi,0,&pidef);
115: mkpf("@e",0,0,0,(int (*)())p_e,const_e,0,&edef);
116:
117: mkpf("log",0,1,uarg,(int (*)())p_log,log,0,&logdef);
118: mkpf("exp",0,1,uarg,(int (*)())p_exp,exp,0,&expdef);
119: mkpf("pow",0,2,darg,(int (*)())p_pow,pow,(int (*)())simplify_pow,&powdef);
120:
121: mkpf("sin",0,1,uarg,(int (*)())p_sin,sin,0,&sindef);
122: mkpf("cos",0,1,uarg,(int (*)())p_cos,cos,0,&cosdef);
123: mkpf("tan",0,1,uarg,(int (*)())p_tan,tan,0,&tandef);
124: mkpf("asin",0,1,uarg,(int (*)())p_asin,asin,0,&asindef);
125: mkpf("acos",0,1,uarg,(int (*)())p_acos,acos,0,&acosdef);
126: mkpf("atan",0,1,uarg,(int (*)())p_atan,atan,0,&atandef);
127:
128: mkpf("sinh",0,1,uarg,(int (*)())p_sinh,sinh,0,&sinhdef);
129: mkpf("cosh",0,1,uarg,(int (*)())p_cosh,cosh,0,&coshdef);
130: mkpf("tanh",0,1,uarg,(int (*)())p_tanh,tanh,0,&tanhdef);
131: #if !defined(VISUAL)
132: mkpf("asinh",0,1,uarg,(int (*)())p_asinh,asinh,0,&asinhdef);
133: mkpf("acosh",0,1,uarg,(int (*)())p_acosh,acosh,0,&acoshdef);
134: mkpf("atanh",0,1,uarg,(int (*)())p_atanh,atanh,0,&atanhdef);
135: #endif
136: make_exp();
137: make_tri();
138: make_itri();
139: make_hyp();
140: #if !defined(VISUAL)
141: make_ihyp();
142: #endif
143: }
144:
145: void make_exp() {
146: V v;
147: P u,vexp,vlog,vpow;
148: Obj *args;
149:
150: mkpfins(expdef,uarg,&v); MKV(v,vexp);
151: mkpfins(powdef,darg,&v); MKV(v,vpow);
152: mkpfins(logdef,uarg,&v); MKV(v,vlog);
153:
154: /* d/dx(log(x)) = 1/x */
155: OALLOC(logdef->deriv,1); divr(CO,(Obj)ONE,(Obj)x,&logdef->deriv[0]);
156:
157: /* d/dx(exp(x)) = exp(x) */
158: OALLOC(expdef->deriv,1); expdef->deriv[0] = (Obj)vexp;
159:
160: /* d/dy(x^y) = log(x)*x^y */
161: OALLOC(powdef->deriv,2); mulp(CO,vpow,vlog,(P *)&powdef->deriv[1]);
162:
163: /* d/dx(x^y) = y*x^(y-1) */
164: args = (Obj *)ALLOCA(2*sizeof(Obj));
165: args[0] = (Obj)x; subp(CO,y,(P)ONE,(P *)&args[1]);
166: _mkpfins(powdef,args,&v); MKV(v,u);
167: mulr(CO,(Obj)u,(Obj)y,&powdef->deriv[0]);
168: }
169:
170: void make_tri() {
171: V v;
172: P vcos,vsin,vtan,t;
173:
174: mkpfins(cosdef,uarg,&v); MKV(v,vcos);
175: mkpfins(sindef,uarg,&v); MKV(v,vsin);
176: mkpfins(tandef,uarg,&v); MKV(v,vtan);
177:
178: /* d/dx(sin(x)) = cos(x) */
179: OALLOC(sindef->deriv,1); sindef->deriv[0] = (Obj)vcos;
180:
181: /* d/dx(cos(x)) = -sin(x) */
182: OALLOC(cosdef->deriv,1); chsgnp(vsin,(P *)&cosdef->deriv[0]);
183:
184: /* d/dx(tan(x)) = 1+tan(x)^2 */
185: OALLOC(tandef->deriv,1);
186: mulr(CO,(Obj)vtan,(Obj)vtan,(Obj *)&t); addp(CO,(P)ONE,t,(P *)&tandef->deriv[0]);
187: }
188:
189: void make_itri() {
190: P t,xx;
191: Q mtwo;
192: V v;
193: Obj *args;
194:
195: /* d/dx(asin(x)) = (1-x^2)^(-1/2) */
196: OALLOC(asindef->deriv,1);
197: args = (Obj *)ALLOCA(2*sizeof(Obj));
198: mulp(CO,x,x,&xx); subp(CO,(P)ONE,xx,(P *)&args[0]);
199: STOQ(-2,mtwo); divq(ONE,mtwo,(Q *)&args[1]);
200: _mkpfins(powdef,args,&v); MKV(v,t);
201: asindef->deriv[0] = (Obj)t;
202:
203: /* d/dx(acos(x)) = -(1-x^2)^(-1/2) */
204: OALLOC(acosdef->deriv,1); chsgnp((P)asindef->deriv[0],(P *)&acosdef->deriv[0]);
205:
206: /* d/dx(atan(x)) = 1/(x^2+1) */
207: OALLOC(atandef->deriv,1);
208: addp(CO,(P)ONE,xx,&t); divr(CO,(Obj)ONE,(Obj)t,&atandef->deriv[0]);
209: }
210:
211: void make_hyp() {
212: V v;
213: P vcosh,vsinh,vtanh,t;
214:
215: mkpfins(coshdef,uarg,&v); MKV(v,vcosh);
216: mkpfins(sinhdef,uarg,&v); MKV(v,vsinh);
217: mkpfins(tanhdef,uarg,&v); MKV(v,vtanh);
218:
219: /* d/dx(sinh(x)) = cosh(x) */
220: OALLOC(sinhdef->deriv,1); sinhdef->deriv[0] = (Obj)vcosh;
221:
222: /* d/dx(cosh(x)) = sinh(x) */
223: OALLOC(coshdef->deriv,1); coshdef->deriv[0] = (Obj)vsinh;
224:
225: /* d/dx(tanh(x)) = 1-tanh(x)^2 */
226: OALLOC(tanhdef->deriv,1);
227: mulr(CO,(Obj)vtanh,(Obj)vtanh,(Obj *)&t); subp(CO,(P)ONE,t,(P *)&tanhdef->deriv[0]);
228: }
229:
230: void make_ihyp() {
231: P t,xx;
232: Q mtwo;
233: V v;
234: Obj *args;
235:
236: /* d/dx(asinh(x)) = (1+x^2)^(-1/2) */
237: OALLOC(asinhdef->deriv,1);
238: args = (Obj *)ALLOCA(2*sizeof(Obj));
239: mulp(CO,x,x,&xx); addp(CO,(P)ONE,xx,(P *)&args[0]);
240: STOQ(-2,mtwo); divq(ONE,mtwo,(Q *)&args[1]);
241: _mkpfins(powdef,args,&v); MKV(v,t);
242: asinhdef->deriv[0] = (Obj)t;
243:
244: /* d/dx(acosh(x)) = (x^2-1)^(-1/2) */
245: OALLOC(acoshdef->deriv,1);
246: subp(CO,xx,(P)ONE,(P *)&args[0]);
247: _mkpfins(powdef,args,&v); MKV(v,t);
248: acoshdef->deriv[0] = (Obj)t;
249:
250: /* d/dx(atanh(x)) = 1/(1-x^2) */
251: OALLOC(atanhdef->deriv,1);
252: subp(CO,(P)ONE,xx,&t); divr(CO,(Obj)ONE,(Obj)t,&atanhdef->deriv[0]);
253: }
254:
255: void mkpow(vl,a,e,r)
256: VL vl;
257: Obj a;
258: Obj e;
259: Obj *r;
260: {
261: PFINS ins;
262: PFAD ad;
263:
264: ins = (PFINS)CALLOC(1,sizeof(PF)+2*sizeof(struct oPFAD));
265: ins->pf = powdef; ad = ins->ad;
266: ad[0].d = 0; ad[0].arg = a; ad[1].d = 0; ad[1].arg = e;
267: simplify_ins(ins,r);
268: }
269:
270: void simplify_pow(ins,rp)
271: PFINS ins;
272: Obj *rp;
273: {
274: PF pf;
275: PFAD ad;
276: Obj a0,a1;
277: V v;
278: P t;
279:
280: pf = ins->pf; ad = ins->ad; a0 = ad[0].arg; a1 = ad[1].arg;
281: if ( !a1 )
282: *rp = (Obj)ONE;
283: else if ( !a0 )
284: *rp = 0;
285: else if ( NUM(a1) && INT(a1) )
286: arf_pwr(CO,a0,a1,rp);
287: else {
288: instov(ins,&v); MKV(v,t); *rp = (Obj)t;
289: }
290: }
291:
292: #define ISPFINS(p)\
293: (p)&&(ID(p) == O_P)&&((int)VR((P)p)->attr!=V_PF)&&\
294: UNIQ(DEG(DC((P)p)))&&UNIQ(COEF(DC((P)p)))
295:
296: void Pfunctor(arg,rp)
297: NODE arg;
298: P *rp;
299: {
300: P p;
301: FUNC t;
302: PF pf;
303: PFINS ins;
304:
305: p = (P)ARG0(arg);
306: if ( !ISPFINS(p) )
307: *rp = 0;
308: else {
309: ins = (PFINS)VR(p)->priv; pf = ins->pf;
310: t = (FUNC)MALLOC(sizeof(struct oFUNC));
311: t->name = pf->name; t->id = A_PURE; t->argc = pf->argc;
312: t->f.puref = pf;
313: makesrvar(t,rp);
314: }
315: }
316:
317: void Pargs(arg,rp)
318: NODE arg;
319: LIST *rp;
320: {
321: P p;
322: PF pf;
323: PFAD ad;
324: PFINS ins;
325: NODE n,n0;
326: int i;
327:
328: p = (P)ARG0(arg);
329: if ( !ISPFINS(p) )
330: *rp = 0;
331: else {
332: ins = (PFINS)VR(p)->priv; ad = ins->ad; pf = ins->pf;
333: for ( i = 0, n0 = 0; i < pf->argc; i++ ) {
334: NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
335: }
336: if ( n0 )
337: NEXT(n) = 0;
338: MKLIST(*rp,n0);
339: }
340: }
341:
342: void Pfunargs(arg,rp)
343: NODE arg;
344: LIST *rp;
345: {
346: P p;
347: P f;
348: FUNC t;
349: PF pf;
350: PFINS ins;
351: PFAD ad;
352: NODE n,n0;
353: int i;
354:
355: p = (P)ARG0(arg);
356: if ( !ISPFINS(p) )
357: *rp = 0;
358: else {
359: ins = (PFINS)VR(p)->priv; ad = ins->ad; pf = ins->pf;
360: t = (FUNC)MALLOC(sizeof(struct oFUNC));
361: t->name = pf->name; t->id = A_PURE; t->argc = pf->argc;
362: t->f.puref = pf;
363: makesrvar(t,&f);
364: n0 = 0; NEXTNODE(n0,n); BDY(n) = (pointer)f;
365: for ( i = 0; i < pf->argc; i++ ) {
366: NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
367: }
368: NEXT(n) = 0;
369: MKLIST(*rp,n0);
370: }
371: }
372:
373: void Pvtype(arg,rp)
374: NODE arg;
375: Q *rp;
376: {
377: P p;
378:
379: p = (P)ARG0(arg);
380: if ( !p || ID(p) != O_P )
381: *rp = 0;
382: else
383: STOQ((int)VR(p)->attr,*rp);
384: }
385:
386: extern FUNC registered_handler;
387:
388: void Pregister_handler(arg,rp)
389: NODE arg;
390: Q *rp;
391: {
392: P p;
393: V v;
394: FUNC func;
395:
396: p = (P)ARG0(arg);
397: if ( !p )
398: registered_handler = 0;
399: else if ( OID(p) != 2 )
400: error("register_hanlder : invalid argument");
401: v = VR(p);
402: if ( (int)v->attr != V_SR )
403: error("register_hanlder : no such function");
404: else {
405: func = (FUNC)v->priv;
406: if ( func->argc )
407: error("register_hanlder : the function must be with no argument");
408: else {
409: registered_handler = func;
410: *rp = ONE;
411: }
412: }
413: }
414:
415: void Pcall(arg,rp)
416: NODE arg;
417: Obj *rp;
418: {
419: P p;
420: V v;
421:
422: p = (P)ARG0(arg);
423: if ( !p || OID(p) != 2 )
424: error("call : invalid argument");
425: v = VR(p);
426: if ( (int)v->attr != V_SR )
427: error("call : no such function");
428:
429: else
430: *rp = (Obj)bevalf((FUNC)v->priv,BDY((LIST)ARG1(arg)));
431: }
432:
433: void Pdeval(arg,rp)
434: NODE arg;
435: Obj *rp;
436: {
437: asir_assert(ARG0(arg),O_R,"deval");
438: devalr(CO,(Obj)ARG0(arg),rp);
439: }
440:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pf.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin