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