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