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