Annotation of OpenXM_contrib2/asir2000/parse/puref.c, Revision 1.15
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.15 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/parse/puref.c,v 1.14 2018/03/28 05:27:22 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
1.13 noro 53: void instoobj(PFINS ins,Obj *rp);
54:
1.1 noro 55: NODE pflist;
56:
1.4 noro 57: void searchpf(char *name,FUNC *fp)
1.1 noro 58: {
1.15 ! noro 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;
1.1 noro 72: }
73:
1.4 noro 74: void searchc(char *name,FUNC *fp)
1.1 noro 75: {
1.15 ! noro 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;
1.1 noro 90: }
91:
1.4 noro 92: void mkpf(char *name,Obj body,int argc,V *args,
1.15 ! noro 93: int (*parif)(),double (*libmf)(), int (*simp)(),PF *pfp)
1.1 noro 94: {
1.15 ! noro 95: PF pf;
! 96: NODE node;
1.1 noro 97:
1.15 ! noro 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;
1.1 noro 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:
1.4 noro 118: void mkpfins(PF pf,V *args,V *vp)
1.1 noro 119: {
1.15 ! noro 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);
1.1 noro 135: }
136:
137: /* the same as above. Argements are given as an array of Obj */
138:
1.4 noro 139: void _mkpfins(PF pf,Obj *args,V *vp)
1.1 noro 140: {
1.15 ! noro 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);
1.1 noro 155: }
156:
157: /* the same as above. darray is also given */
158:
1.4 noro 159: void _mkpfins_with_darray(PF pf,Obj *args,int *darray,V *vp)
1.1 noro 160: {
1.15 ! noro 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);
1.1 noro 175: }
176:
1.4 noro 177: void appendpfins(V v,V *vp)
1.1 noro 178: {
1.15 ! noro 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;
1.1 noro 197: }
198:
1.4 noro 199: void duppfins(V v,V *vp)
1.1 noro 200: {
1.15 ! noro 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;
1.1 noro 210: }
211:
1.4 noro 212: void derivvar(VL vl,V pf,V v,Obj *a)
1.1 noro 213: {
1.15 ! noro 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: }
1.1 noro 240: }
241:
1.4 noro 242: void derivr(VL vl,Obj a,V v,Obj *b)
1.1 noro 243: {
1.15 ! noro 244: VL tvl,svl;
! 245: Obj r,s,t,u,nm,dn,dnm,ddn,m;
1.1 noro 246:
1.15 ! noro 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: }
1.8 noro 287: }
288:
1.9 noro 289: void simple_derivr(VL vl,Obj a,V v,Obj *b)
290: {
1.15 ! noro 291: Obj r,s,t,u,nm,dn;
1.9 noro 292:
1.15 ! noro 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: }
1.9 noro 309: }
310:
1.8 noro 311: int obj_is_dependent(Obj a,V v)
312: {
1.15 ! noro 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");
1.8 noro 319: }
320:
321: int poly_is_dependent(P p,V v)
322: {
1.15 ! noro 323: DCP dc;
1.8 noro 324:
1.15 ! noro 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: }
1.1 noro 332: }
333:
1.7 noro 334: void gen_pwrr(VL vl,Obj a,Obj d,Obj *r)
335: {
1.15 ! noro 336: if ( INT(d) )
! 337: pwrr(vl,a,d,r);
! 338: else
! 339: mkpow(vl,a,d,r);
1.7 noro 340: }
341:
1.4 noro 342: void substr(VL vl,int partial,Obj a,V v,Obj b,Obj *c)
1.1 noro 343: {
1.15 ! noro 344: Obj nm,dn,t;
1.1 noro 345:
1.15 ! noro 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: }
1.1 noro 363: }
364:
1.4 noro 365: void substpr(VL vl,int partial,Obj p,V v0,Obj p0,Obj *pr)
1.1 noro 366: {
1.15 ! noro 367: P x;
! 368: Obj t,m,c,s,a;
! 369: DCP dc;
! 370: Q 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: subq(d,DEG(dc),(Q *)&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: }
1.1 noro 426: }
427:
1.4 noro 428: void evalr(VL vl,Obj a,int prec,Obj *c)
1.1 noro 429: {
1.15 ! noro 430: Obj nm,dn;
1.1 noro 431:
1.15 ! noro 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: }
1.1 noro 448: }
449:
1.4 noro 450: void evalp(VL vl,P p,int prec,P *pr)
1.1 noro 451: {
1.15 ! noro 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: }
1.1 noro 476: }
477:
1.4 noro 478: void evalv(VL vl,V v,int prec,Obj *rp)
1.1 noro 479: {
1.15 ! noro 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: }
1.1 noro 497: }
498:
1.4 noro 499: void evalins(PFINS ins,int prec,Obj *rp)
1.1 noro 500: {
1.15 ! noro 501: PF pf;
! 502: PFINS tins;
! 503: PFAD ad,tad;
! 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: 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 ) {
! 525: NEXTNODE(n0,n); STOQ(prec,q); BDY(n) = (pointer)q;
! 526: }
! 527: if ( n0 )
! 528: NEXT(n) = 0;
! 529: (*pf->pari)(n0,rp);
! 530: }
1.1 noro 531: }
532:
1.4 noro 533: void devalr(VL vl,Obj a,Obj *c)
1.1 noro 534: {
1.15 ! noro 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(re,im,&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: }
1.1 noro 570: }
571:
1.4 noro 572: void devalp(VL vl,P p,P *pr)
1.1 noro 573: {
1.15 ! noro 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: }
1.1 noro 613: }
614:
1.4 noro 615: void devalv(VL vl,V v,Obj *rp)
1.1 noro 616: {
1.15 ! noro 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: }
1.1 noro 634: }
635:
1.4 noro 636: void devalins(PFINS ins,Obj *rp)
1.1 noro 637: {
1.15 ! noro 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: }
1.1 noro 679: }
680:
1.13 noro 681: extern int evalef,bigfloat;
1.10 noro 682:
683: void simplify_elemfunc_ins(PFINS ins,Obj *rp)
684: {
1.13 noro 685: if ( evalef ) {
686: if ( bigfloat ) evalins(ins,0,rp);
687: else devalins(ins,rp);
688: } else instoobj(ins,rp);
1.10 noro 689: }
690:
1.14 noro 691: void simplify_factorial_ins(PFINS ins,Obj *rp)
692: {
693: PFAD ad;
694: Obj a;
695: Q q;
696:
697: ad = ins->ad;
698: a = ad[0].arg;
699: if ( !ad[0].d && INT(a) && ( !a || (PL(NM((Q)a)) == 1 && SGN((Q)a) > 0) ) ) {
700: factorial(QTOS((Q)a),&q);
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 || SGN((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,rp);
727: } else simplify_elemfunc_ins(ins,rp);
728: }
729:
1.4 noro 730: void simplify_ins(PFINS ins,Obj *rp)
1.1 noro 731: {
1.15 ! noro 732: V v;
! 733: P t;
1.1 noro 734:
1.15 ! noro 735: if ( ins->pf->simplify )
! 736: (*ins->pf->simplify)(ins,rp);
! 737: else {
! 738: instoobj(ins,rp);
! 739: }
1.1 noro 740: }
741:
1.13 noro 742: void instoobj(PFINS ins,Obj *rp)
1.1 noro 743: {
1.15 ! noro 744: V v,newv;
! 745: P t;
1.1 noro 746:
1.15 ! noro 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;
1.1 noro 752: }
753:
1.4 noro 754: void substfr(VL vl,Obj a,PF u,PF f,Obj *c)
1.1 noro 755: {
1.15 ! noro 756: Obj nm,dn;
1.1 noro 757:
1.15 ! noro 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: }
1.1 noro 774: }
775:
1.4 noro 776: void substfp(VL vl,Obj p,PF u,PF f,Obj *pr)
1.1 noro 777: {
1.15 ! noro 778: V v;
! 779: DCP dc;
! 780: Obj a,c,m,s,t,p0;
! 781: Q 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 ( (int)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: subq(d,DEG(dc),(Q *)&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: }
1.1 noro 817: }
818:
1.4 noro 819: void substfv(VL vl,V v,PF u,PF f,Obj *c)
1.1 noro 820: {
1.15 ! noro 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: }
1.1 noro 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 / asir2000 / parse