Annotation of OpenXM_contrib2/asir2018/parse/puref.c, Revision 1.1
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: *
! 48: * $OpenXM$
! 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 ) {
! 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: }
! 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) ) {
! 700: factorialz(QTOS((Z)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 || 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