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