![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to puref.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / parse|
Return to puref.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / parse |
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/parse/puref.c,v 1.1.1.1 1999/12/03 07:39:12 noro Exp $
! 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
53: NODE pflist;
54:
55: void searchpf(name,fp)
56: char *name;
57: FUNC *fp;
58: {
59: NODE node;
60: PF pf;
61: FUNC t;
62:
63: for ( node = pflist; node; node = NEXT(node) )
64: if ( !strcmp(name,((PF)node->body)->name) ) {
65: pf = (PF)node->body;
66: *fp = t = (FUNC)MALLOC(sizeof(struct oFUNC));
67: t->name = name; t->id = A_PURE; t->argc = pf->argc;
68: t->f.puref = pf;
69: return;
70: }
71: *fp = 0;
72: }
73:
74: void searchc(name,fp)
75: char *name;
76: FUNC *fp;
77: {
78: NODE node;
79: PF pf;
80: FUNC t;
81:
82: for ( node = pflist; node; node = NEXT(node) )
83: if ( !strcmp(name,((PF)node->body)->name)
84: && !((PF)node->body)->argc ) {
85: pf = (PF)node->body;
86: *fp = t = (FUNC)MALLOC(sizeof(struct oFUNC));
87: t->name = name; t->id = A_PURE; t->argc = pf->argc;
88: t->f.puref = pf;
89: return;
90: }
91: *fp = 0;
92: }
93:
94: void mkpf(name,body,argc,args,parif,libmf,simp,pfp)
95: char *name;
96: Obj body;
97: int argc;
98: V *args;
99: int (*parif)(),(*simp)();
100: double (*libmf)();
101: PF *pfp;
102: {
103: PF pf;
104: NODE node;
105:
106: NEWPF(pf); pf->name = name; pf->body = body;
107: pf->argc = argc; pf->args = args; pf->pari = parif; pf->simplify = simp;
108: pf->libm = libmf;
109: for ( node = pflist; node; node = NEXT(node) )
110: if ( !strcmp(((PF)BDY(node))->name,name) )
111: break;
112: if ( !node ) {
113: NEWNODE(node); NEXT(node) = pflist; pflist = node;
114: /* fprintf(stderr,"%s() defined.\n",name); */
115: } else
116: fprintf(stderr,"%s() redefined.\n",name);
117: BDY(node) = (pointer)pf; *pfp = pf;
118: }
119:
120: /*
121: create an instance of a pure function. args are given
122: as an array of V. So we have to create a P object for
123: each arg.
124: */
125:
126: void mkpfins(pf,args,vp)
127: PF pf;
128: V *args;
129: V *vp;
130: {
131: V v;
132: PFINS ins;
133: PFAD ad;
134: int i;
135: P t;
136:
137: NEWV(v); NAME(v) = 0; v->attr = (pointer)V_PF;
138: ins = (PFINS)MALLOC(sizeof(PF)+pf->argc*sizeof(struct oPFAD));
139: bzero((char *)ins,(int)(sizeof(PF)+pf->argc*sizeof(struct oPFAD)));
140: ins->pf = pf;
141: v->priv = (pointer)ins;
142: for ( i = 0, ad = ins->ad; i < pf->argc; i++ ) {
143: ad[i].d = 0; MKV(args[i],t); ad[i].arg = (Obj)t;
144: }
145: appendpfins(v,vp);
146: }
147:
148: /* the same as above. Argements are given as an array of Obj */
149:
150: void _mkpfins(pf,args,vp)
151: PF pf;
152: Obj *args;
153: V *vp;
154: {
155: V v;
156: PFINS ins;
157: PFAD ad;
158: int i;
159:
160: NEWV(v); NAME(v) = 0; v->attr = (pointer)V_PF;
161: ins = (PFINS)MALLOC(sizeof(PF)+pf->argc*sizeof(struct oPFAD));
162: bzero((char *)ins,(int)(sizeof(PF)+pf->argc*sizeof(struct oPFAD)));
163: ins->pf = pf;
164: v->priv = (pointer)ins;
165: for ( i = 0, ad = ins->ad; i < pf->argc; i++ ) {
166: ad[i].d = 0; ad[i].arg = args[i];
167: }
168: appendpfins(v,vp);
169: }
170:
171: /* the same as above. darray is also given */
172:
173: void _mkpfins_with_darray(pf,args,darray,vp)
174: PF pf;
175: Obj *args;
176: int *darray;
177: V *vp;
178: {
179: V v;
180: PFINS ins;
181: PFAD ad;
182: int i;
183:
184: NEWV(v); NAME(v) = 0; v->attr = (pointer)V_PF;
185: ins = (PFINS)MALLOC(sizeof(PF)+pf->argc*sizeof(struct oPFAD));
186: bzero((char *)ins,(int)(sizeof(PF)+pf->argc*sizeof(struct oPFAD)));
187: ins->pf = pf;
188: v->priv = (pointer)ins;
189: for ( i = 0, ad = ins->ad; i < pf->argc; i++ ) {
190: ad[i].d = darray[i]; ad[i].arg = args[i];
191: }
192: appendpfins(v,vp);
193: }
194:
195: void appendpfins(v,vp)
196: V v;
197: V *vp;
198: {
199: PF fdef;
200: PFAD ad,tad;
201: NODE node;
202: int i;
203:
204: fdef = ((PFINS)v->priv)->pf; ad = ((PFINS)v->priv)->ad;
205: for ( node = fdef->ins; node; node = NEXT(node) ) {
206: for ( i = 0, tad = ((PFINS)((V)node->body)->priv)->ad;
207: i < fdef->argc; i++ )
208: if ( (ad[i].d != tad[i].d) || compr(CO,ad[i].arg,tad[i].arg) )
209: break;
210: if ( i == fdef->argc ) {
211: *vp = (V)node->body;
212: return;
213: }
214: }
215: NEWNODE(node); node->body = (pointer)v; NEXT(node) = fdef->ins;
216: fdef->ins = node; appendvar(CO,v); *vp = v;
217: }
218:
219: void duppfins(v,vp)
220: V v;
221: V *vp;
222: {
223: V tv;
224: PFINS tins;
225: int size;
226:
227: NEWV(tv); tv->name = v->name; tv->attr = v->attr;
228: size = sizeof(PF)+((PFINS)v->priv)->pf->argc*sizeof(struct oPFAD);
229: tins = (PFINS)MALLOC(size); bcopy((char *)v->priv,(char *)tins,size);
230: tv->priv = (pointer)tins;
231: *vp = tv;
232: }
233:
234: void derivvar(vl,pf,v,a)
235: VL vl;
236: V pf,v;
237: Obj *a;
238: {
239: Obj t,s,u,w,u1;
240: P p;
241: V tv,sv;
242: PF fdef;
243: PFAD ad;
244: int i,j;
245:
246: fdef = ((PFINS)pf->priv)->pf; ad = ((PFINS)pf->priv)->ad;
247: if ( fdef->deriv ) {
248: for ( t = 0, i = 0; i < fdef->argc; i++ ) {
249: derivr(vl,ad[i].arg,v,&s);
250: for ( j = 0, u = fdef->deriv[i]; j < fdef->argc; j++ ) {
251: substr(vl,0,u,fdef->args[j],ad[j].arg,&u1); u = u1;
252: }
253: mulr(vl,s,u,&w); addr(vl,t,w,&s); t = s;
254: }
255: *a = t;
256: } else {
257: for ( t = 0, i = 0; i < fdef->argc; i++ ) {
258: derivr(vl,ad[i].arg,v,&s);
259: duppfins(pf,&tv); (((PFINS)tv->priv)->ad)[i].d++;
260: appendpfins(tv,&sv);
261: MKV(sv,p); mulr(vl,s,(Obj)p,&w); addr(vl,t,w,&s); t = s;
262: }
263: *a = t;
264: }
265: }
266:
267: void derivr(vl,a,v,b)
268: VL vl;
269: V v;
270: Obj a,*b;
271: {
272: VL tvl,svl;
273: Obj r,s,t,u,nm,dn,dnm,ddn,m;
274:
275: if ( !a )
276: *b = 0;
277: else
278: switch ( OID(a) ) {
279: case O_N:
280: *b = 0; break;
281: case O_P:
282: clctvr(vl,a,&tvl);
283: for ( dnm = 0, svl = tvl; svl; svl = NEXT(svl) ) {
284: if ( svl->v == v ) {
285: pderivr(vl,a,v,&s); addr(vl,s,dnm,&u); dnm = u;
286: } else if ( (vid)svl->v->attr == V_PF ) {
287: pderivr(vl,a,svl->v,&s); derivvar(vl,svl->v,v,&r);
288: mulr(vl,s,r,&u); addr(vl,u,dnm,&s); dnm = s;
289: }
290: }
291: *b = (Obj)dnm; break;
292: case O_R:
293: clctvr(vl,a,&tvl);
294: nm = (Obj)NM((R)a); dn = (Obj)DN((R)a);
295: for ( dnm = ddn = 0, svl = tvl; svl; svl = NEXT(svl) ) {
296: if ( svl->v == v ) {
297: pderivr(vl,nm,v,&s); addr(vl,s,dnm,&u); dnm = u;
298: pderivr(vl,dn,v,&s); addr(vl,s,ddn,&u); ddn = u;
299: } else if ( (vid)svl->v->attr == V_PF ) {
300: pderivr(vl,nm,svl->v,&s); derivvar(vl,svl->v,v,&r);
301: mulr(vl,s,r,&u); addr(vl,u,dnm,&s); dnm = s;
302: pderivr(vl,dn,svl->v,&s); derivvar(vl,svl->v,v,&r);
303: mulr(vl,s,r,&u); addr(vl,u,ddn,&s); ddn = s;
304: }
305: }
306: mulr(vl,dnm,dn,&t); mulr(vl,ddn,nm,&s);
307: subr(vl,t,s,&u); reductr(vl,u,&t);
308: if ( !t )
309: *b = 0;
310: else {
311: mulp(vl,(P)dn,(P)dn,(P *)&m); divr(vl,t,m,b);
312: }
313: break;
314: }
315: }
316:
317: void substr(vl,partial,a,v,b,c)
318: VL vl;
319: int partial;
320: Obj a;
321: V v;
322: Obj b;
323: Obj *c;
324: {
325: Obj nm,dn,t;
326:
327: if ( !a )
328: *c = 0;
329: else {
330: switch ( OID(a) ) {
331: case O_N:
332: *c = a; break;
333: case O_P:
334: substpr(vl,partial,a,v,b,c); break;
335: case O_R:
336: substpr(vl,partial,(Obj)NM((R)a),v,b,&nm);
337: substpr(vl,partial,(Obj)DN((R)a),v,b,&dn);
338: divr(vl,nm,dn,&t); reductr(vl,t,c);
339: break;
340: default:
341: *c = 0; break;
342: }
343: }
344: }
345:
346: void substpr(vl,partial,p,v0,p0,pr)
347: VL vl;
348: int partial;
349: V v0;
350: Obj p;
351: Obj p0;
352: Obj *pr;
353: {
354: P x;
355: Obj t,m,c,s,a;
356: DCP dc;
357: Q d;
358: V v;
359: PF pf;
360: PFAD ad,tad;
361: PFINS tins;
362: int i;
363:
364: if ( !p )
365: *pr = 0;
366: else if ( NUM(p) )
367: *pr = (Obj)p;
368: else if ( (v = VR((P)p)) != v0 ) {
369: if ( !partial && ((vid)v->attr == V_PF) ) {
370: ad = ((PFINS)v->priv)->ad; pf = ((PFINS)v->priv)->pf;
371: tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
372: tins->pf = pf;
373: for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) {
374: tad[i].d = ad[i].d;
375: substr(vl,partial,ad[i].arg,v0,p0,&tad[i].arg);
376: }
377: simplify_ins(tins,(Obj *)&x);
378: } else
379: MKV(VR((P)p),x);
380: for ( c = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
381: substpr(vl,partial,(Obj)COEF(dc),v0,p0,&t);
382: if ( DEG(dc) ) {
383: pwrp(vl,x,DEG(dc),(P *)&s); mulr(vl,s,t,&m);
384: addr(vl,m,c,&a); c = a;
385: } else {
386: addr(vl,t,c,&a); c = a;
387: }
388: }
389: *pr = c;
390: } else {
391: dc = DC((P)p);
392: if ( !partial )
393: substpr(vl,partial,(Obj)COEF(dc),v0,p0,&c);
394: else
395: c = (Obj)COEF(dc);
396: for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc = NEXT(dc) ) {
397: subq(d,DEG(dc),(Q *)&t); pwrr(vl,p0,t,&s); mulr(vl,s,c,&m);
398: if ( !partial )
399: substpr(vl,partial,(Obj)COEF(dc),v0,p0,&t);
400: else
401: t = (Obj)COEF(dc);
402: addr(vl,m,t,&c);
403: }
404: if ( d ) {
405: pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m);
406: c = m;
407: }
408: *pr = c;
409: }
410: }
411:
412: void evalr(vl,a,prec,c)
413: VL vl;
414: Obj a;
415: int prec;
416: Obj *c;
417: {
418: Obj nm,dn;
419:
420: if ( !a )
421: *c = 0;
422: else {
423: switch ( OID(a) ) {
424: case O_N:
425: *c = a; break;
426: case O_P:
427: evalp(vl,(P)a,prec,(P *)c); break;
428: case O_R:
429: evalp(vl,NM((R)a),prec,(P *)&nm); evalp(vl,DN((R)a),prec,(P *)&dn);
430: divr(vl,nm,dn,c);
431: break;
432: default:
433: error("evalr : not implemented"); break;
434: }
435: }
436: }
437:
438: void evalp(vl,p,prec,pr)
439: VL vl;
440: P p;
441: int prec;
442: P *pr;
443: {
444: P t;
445: DCP dc,dcr0,dcr;
446: Obj u;
447:
448: if ( !p || NUM(p) )
449: *pr = p;
450: else {
451: for ( dcr0 = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
452: evalp(vl,COEF(dc),prec,&t);
453: if ( t ) {
454: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
455: }
456: }
457: if ( !dcr0 ) {
458: *pr = 0; return;
459: } else {
460: NEXT(dcr) = 0; MKP(VR(p),dcr0,t);
461: }
462: if ( NUM(t) || (VR(t) != VR(p)) || ((vid)VR(p)->attr != V_PF) ) {
463: *pr = t; return;
464: } else {
465: evalv(vl,VR(p),prec,&u); substr(vl,1,(Obj)t,VR(p),u,(Obj *)pr);
466: }
467: }
468: }
469:
470: void evalv(vl,v,prec,rp)
471: VL vl;
472: V v;
473: int prec;
474: Obj *rp;
475: {
476: PFINS ins,tins;
477: PFAD ad,tad;
478: PF pf;
479: P t;
480: int i;
481:
482: if ( (vid)v->attr != V_PF ) {
483: MKV(v,t); *rp = (Obj)t;
484: } else {
485: ins = (PFINS)v->priv; ad = ins->ad; pf = ins->pf;
486: tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
487: tins->pf = pf;
488: for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) {
489: tad[i].d = ad[i].d; evalr(vl,ad[i].arg,prec,&tad[i].arg);
490: }
491: evalins(tins,prec,rp);
492: }
493: }
494:
495: void evalins(ins,prec,rp)
496: PFINS ins;
497: int prec;
498: Obj *rp;
499: {
500: PF pf;
501: PFAD ad;
502: int i;
503: Q q;
504: V v;
505: P x;
506: NODE n0,n;
507:
508: pf = ins->pf; ad = ins->ad;
509: for ( i = 0; i < pf->argc; i++ )
510: if ( ad[i].d || (ad[i].arg && !NUM(ad[i].arg)) )
511: break;
512: if ( (i != pf->argc) || !pf->pari ) {
513: instov(ins,&v); MKV(v,x); *rp = (Obj)x;
514: } else {
515: for ( n0 = 0, i = 0; i < pf->argc; i++ ) {
516: NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
517: }
518: if ( prec ) {
519: NEXTNODE(n0,n); STOQ(prec,q); BDY(n) = (pointer)q;
520: }
521: if ( n0 )
522: NEXT(n) = 0;
523: (*pf->pari)(n0,rp);
524: }
525: }
526:
527: void devalins(PFINS,Obj *);
528: void devalv(VL,V,Obj *);
529: void devalp(VL,P,P *);
530: void devalr(VL,Obj,Obj *);
531:
532: void devalr(vl,a,c)
533: VL vl;
534: Obj a;
535: Obj *c;
536: {
537: Obj nm,dn;
538: double d;
539: Real r;
540:
541: if ( !a )
542: *c = 0;
543: else {
544: switch ( OID(a) ) {
545: case O_N:
546: d = ToReal(a);
547: MKReal(d,r);
548: *c = (Obj)r;
549: break;
550: case O_P:
551: devalp(vl,(P)a,(P *)c); break;
552: case O_R:
553: devalp(vl,NM((R)a),(P *)&nm);
554: devalp(vl,DN((R)a),(P *)&dn);
555: divr(vl,nm,dn,c);
556: break;
557: default:
558: error("devalr : not implemented"); break;
559: }
560: }
561: }
562:
563: void devalp(vl,p,pr)
564: VL vl;
565: P p;
566: P *pr;
567: {
568: P t;
569: DCP dc,dcr0,dcr;
570: Obj u,s;
571: double d;
572: Real r;
573:
574: if ( !p || NUM(p) ) {
575: d = ToReal(p);
576: MKReal(d,r);
577: *pr = (P)r;
578: } else {
579: for ( dcr0 = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
580: devalp(vl,COEF(dc),&t);
581: if ( t ) {
582: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
583: }
584: }
585: if ( !dcr0 )
586: *pr = 0;
587: else {
588: NEXT(dcr) = 0; MKP(VR(p),dcr0,t);
589: if ( NUM(t) ) {
590: d = ToReal((Num)t);
591: MKReal(d,r);
592: *pr = (P)r;
593: } else if ( (VR(t) != VR(p)) || (VR(p)->attr != (pointer)V_PF) )
594: *pr = t;
595: else {
596: devalv(vl,VR(p),&u);
597: substr(vl,1,(Obj)t,VR(p),u,&s);
598: if ( s && NUM(s) ) {
599: d = ToReal((Num)s);
600: MKReal(d,r);
601: *pr = (P)r;
602: } else
603: *pr = (P)s;
604: }
605: }
606: }
607: }
608:
609: void devalv(vl,v,rp)
610: VL vl;
611: V v;
612: Obj *rp;
613: {
614: PFINS ins,tins;
615: PFAD ad,tad;
616: PF pf;
617: P t;
618: Obj s;
619: int i;
620:
621: if ( (vid)v->attr != V_PF ) {
622: MKV(v,t); *rp = (Obj)t;
623: } else {
624: ins = (PFINS)v->priv; ad = ins->ad; pf = ins->pf;
625: tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
626: tins->pf = pf;
627: for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) {
628: tad[i].d = ad[i].d; devalr(vl,ad[i].arg,&tad[i].arg);
629: }
630: devalins(tins,rp);
631: }
632: }
633:
634: void devalins(ins,rp)
635: PFINS ins;
636: Obj *rp;
637: {
638: PF pf;
639: PFAD ad;
640: int i;
641: Real r;
642: double d;
643: Q q;
644: V v;
645: P x;
646:
647: pf = ins->pf; ad = ins->ad;
648: for ( i = 0; i < pf->argc; i++ )
649: if ( ad[i].d || (ad[i].arg && !NUM(ad[i].arg)) )
650: break;
651: if ( (i != pf->argc) || !pf->libm ) {
652: instov(ins,&v); MKV(v,x); *rp = (Obj)x;
653: } else {
654: switch ( pf->argc ) {
655: case 0:
656: d = (*pf->libm)(); break;
657: case 1:
658: d = (*pf->libm)(ToReal(ad[0].arg)); break;
659: case 2:
660: d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg)); break;
661: case 3:
662: d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg),
663: ToReal(ad[2].arg)); break;
664: case 4:
665: d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg),
666: ToReal(ad[2].arg),ToReal(ad[3].arg)); break;
667: default:
668: error("devalv : not supported");
669: }
670: MKReal(d,r); *rp = (Obj)r;
671: }
672: }
673:
674: void simplify_ins(ins,rp)
675: PFINS ins;
676: Obj *rp;
677: {
678: V v;
679: P t;
680:
681: if ( ins->pf->simplify )
682: (*ins->pf->simplify)(ins,rp);
683: else {
684: instov(ins,&v); MKV(v,t); *rp = (Obj)t;
685: }
686: }
687:
688: void instov(ins,vp)
689: PFINS ins;
690: V *vp;
691: {
692: V v;
693:
694: NEWV(v); NAME(v) = 0;
695: v->attr = (pointer)V_PF; v->priv = (pointer)ins;
696: appendpfins(v,vp);
697: }
698:
699: void substfr(vl,a,u,f,c)
700: VL vl;
701: Obj a;
702: PF u,f;
703: Obj *c;
704: {
705: Obj nm,dn;
706:
707: if ( !a )
708: *c = 0;
709: else {
710: switch ( OID(a) ) {
711: case O_N:
712: *c = a; break;
713: case O_P:
714: substfp(vl,a,u,f,c); break;
715: case O_R:
716: substfp(vl,(Obj)NM((R)a),u,f,&nm); substfp(vl,(Obj)DN((R)a),u,f,&dn);
717: divr(vl,nm,dn,c);
718: break;
719: default:
720: error("substfr : not implemented"); break;
721: }
722: }
723: }
724:
725: void substfp(vl,p,u,f,pr)
726: VL vl;
727: Obj p;
728: PF u,f;
729: Obj *pr;
730: {
731: V v;
732: DCP dc;
733: Obj a,c,m,s,t,p0;
734: Q d;
735: P x;
736:
737: if ( !p )
738: *pr = 0;
739: else if ( NUM(p) )
740: *pr = (Obj)p;
741: else {
742: v = VR((P)p); dc = DC((P)p);
743: if ( (int)v->attr != V_PF ) {
744: MKV(VR((P)p),x);
745: for ( c = 0; dc; dc = NEXT(dc) ) {
746: substfp(vl,(Obj)COEF(dc),u,f,&t);
747: if ( DEG(dc) ) {
748: pwrp(vl,x,DEG(dc),(P *)&s); mulr(vl,s,t,&m);
749: addr(vl,m,c,&a); c = a;
750: } else {
751: addr(vl,t,c,&a); c = a;
752: }
753: }
754: } else {
755: substfv(vl,v,u,f,&p0);
756: substfp(vl,(Obj)COEF(dc),u,f,&c);
757: for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc = NEXT(dc) ) {
758: subq(d,DEG(dc),(Q *)&t); pwrr(vl,p0,t,&s); mulr(vl,s,c,&m);
759: substfp(vl,(Obj)COEF(dc),u,f,&t); addr(vl,m,t,&c);
760: }
761: if ( d ) {
762: pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m);
763: c = m;
764: }
765: }
766: *pr = c;
767: }
768: }
769:
770: void substfv(vl,v,u,f,c)
771: VL vl;
772: V v;
773: PF u,f;
774: Obj *c;
775: {
776: P t;
777: Obj r,s,w;
778: int i,j;
779: PFINS ins,tins;
780: PFAD ad,tad;
781:
782: ins = (PFINS)v->priv; ad = ins->ad;
783: if ( ins->pf == u ) {
784: if ( u->argc != f->argc )
785: error("substfv : argument mismatch");
786: if ( !f->body ) {
787: mkpfins(f,f->args,&v); MKV(v,t); r = (Obj)t;
788: } else
789: r = f->body;
790: for ( i = 0; i < f->argc; i++ )
791: for ( j = 0; j < ad[i].d; j++ ) {
792: derivr(vl,r,f->args[i],&s); r = s;
793: }
794: for ( i = 0; i < f->argc; i++ ) {
795: substfr(vl,ad[i].arg,u,f,&w);
796: substr(vl,0,r,f->args[i],w,&s); r = s;
797: }
798: *c = r;
799: } else {
800: tins = (PFINS)MALLOC(sizeof(PF)+f->argc*sizeof(struct oPFAD));
801: tins->pf = ins->pf; tad = tins->ad;
802: for ( i = 0; i < f->argc; i++ ) {
803: tad[i].d = ad[i].d; substfr(vl,ad[i].arg,u,f,&tad[i].arg);
804: }
805: instov(tins,&v); MKV(v,t); *c = (Obj)t;
806: }
807: }