[BACK]Return to puref.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / parse

Annotation of OpenXM_contrib2/asir2018/parse/puref.c, Revision 1.4

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.4     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2018/parse/puref.c,v 1.3 2019/11/12 10:53:23 kondoh 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);
1.4     ! noro      325:   else {
1.1       noro      326:     error("obj_is_dependent : not implemented");
1.4     ! noro      327:     return 0;
        !           328:   }
1.1       noro      329: }
                    330:
                    331: int poly_is_dependent(P p,V v)
                    332: {
                    333:   DCP dc;
                    334:
                    335:   if ( !p || OID(p) <= O_N ) return 0;
                    336:   else if ( v == VR(p) ) return 1;
                    337:   else {
                    338:     for ( dc = DC(p); dc; dc = NEXT(dc) )
                    339:       if ( poly_is_dependent(COEF(dc),v) ) return 1;
                    340:     return 0;
                    341:   }
                    342: }
                    343:
                    344: void gen_pwrr(VL vl,Obj a,Obj d,Obj *r)
                    345: {
                    346:   if ( INT(d) )
                    347:     pwrr(vl,a,d,r);
                    348:   else
                    349:     mkpow(vl,a,d,r);
                    350: }
                    351:
                    352: void substr(VL vl,int partial,Obj a,V v,Obj b,Obj *c)
                    353: {
                    354:   Obj nm,dn,t;
                    355:
                    356:   if ( !a )
                    357:     *c = 0;
                    358:   else {
                    359:     switch ( OID(a) ) {
                    360:       case O_N:
                    361:         *c = a; break;
                    362:       case O_P:
                    363:         substpr(vl,partial,a,v,b,c); break;
                    364:       case O_R:
                    365:         substpr(vl,partial,(Obj)NM((R)a),v,b,&nm);
                    366:         substpr(vl,partial,(Obj)DN((R)a),v,b,&dn);
                    367:         divr(vl,nm,dn,&t); reductr(vl,t,c);
                    368:         break;
                    369:       default:
                    370:         *c = 0; break;
                    371:     }
                    372:   }
                    373: }
                    374:
                    375: void substpr(VL vl,int partial,Obj p,V v0,Obj p0,Obj *pr)
                    376: {
                    377:   P x;
                    378:   Obj t,m,c,s,a;
                    379:   DCP dc;
                    380:   Z d;
                    381:   V v;
                    382:   PF pf;
                    383:   PFAD ad,tad;
                    384:   PFINS tins;
                    385:   int i;
                    386:
                    387:   if ( !p )
                    388:     *pr = 0;
                    389:   else if ( NUM(p) )
                    390:     *pr = (Obj)p;
                    391:   else if ( (v = VR((P)p)) != v0 ) {
                    392:     if ( !partial && ((vid)v->attr == V_PF) ) {
                    393:       ad = ((PFINS)v->priv)->ad; pf = ((PFINS)v->priv)->pf;
                    394:       tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
                    395:       tins->pf = pf;
                    396:       for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) {
                    397:         tad[i].d = ad[i].d;
                    398:         substr(vl,partial,ad[i].arg,v0,p0,&tad[i].arg);
                    399:       }
                    400:       simplify_ins(tins,(Obj *)&x);
                    401:     } else
                    402:       MKV(VR((P)p),x);
                    403:     for ( c = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
                    404:       substpr(vl,partial,(Obj)COEF(dc),v0,p0,&t);
                    405:       if ( DEG(dc) ) {
                    406:         gen_pwrr(vl,(Obj)x,(Obj)DEG(dc),&s);
                    407:         mulr(vl,s,t,&m);
                    408:         addr(vl,m,c,&a); c = a;
                    409:       } else {
                    410:         addr(vl,t,c,&a); c = a;
                    411:       }
                    412:     }
                    413:     *pr = c;
                    414:   } else {
                    415:     dc = DC((P)p);
                    416:     if ( !partial )
                    417:       substpr(vl,partial,(Obj)COEF(dc),v0,p0,&c);
                    418:     else
                    419:       c = (Obj)COEF(dc);
                    420:     for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc = NEXT(dc) ) {
                    421:         subz(d,DEG(dc),(Z *)&t);
                    422:         gen_pwrr(vl,p0,t,&s); mulr(vl,s,c,&m);
                    423:         if ( !partial )
                    424:           substpr(vl,partial,(Obj)COEF(dc),v0,p0,&t);
                    425:         else
                    426:           t = (Obj)COEF(dc);
                    427:         addr(vl,m,t,&c);
                    428:     }
                    429:     if ( d ) {
                    430:       gen_pwrr(vl,p0,(Obj)d,&t);
                    431:       mulr(vl,t,c,&m);
                    432:       c = m;
                    433:     }
                    434:     *pr = c;
                    435:   }
                    436: }
                    437:
                    438: void evalr(VL vl,Obj a,int prec,Obj *c)
                    439: {
                    440:   Obj nm,dn;
                    441:
                    442:   if ( !a )
                    443:     *c = 0;
                    444:   else {
                    445:     switch ( OID(a) ) {
                    446:       case O_N:
                    447:         *c = a; break;
                    448:       case O_P:
                    449:         evalp(vl,(P)a,prec,(P *)c); break;
                    450:       case O_R:
                    451:         evalp(vl,NM((R)a),prec,(P *)&nm); evalp(vl,DN((R)a),prec,(P *)&dn);
                    452:         divr(vl,nm,dn,c);
                    453:         break;
                    454:       default:
                    455:         error("evalr : not implemented"); break;
                    456:     }
                    457:   }
                    458: }
                    459:
                    460: void evalp(VL vl,P p,int prec,P *pr)
                    461: {
                    462:   P t;
                    463:   DCP dc,dcr0,dcr;
                    464:   Obj u;
                    465:
                    466:   if ( !p || NUM(p) )
                    467:     *pr = p;
                    468:   else {
                    469:     for ( dcr0 = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
                    470:       evalp(vl,COEF(dc),prec,&t);
                    471:       if ( t ) {
                    472:         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
                    473:       }
                    474:     }
                    475:     if ( !dcr0 ) {
                    476:       *pr = 0; return;
                    477:     } else {
                    478:       NEXT(dcr) = 0; MKP(VR(p),dcr0,t);
                    479:     }
                    480:     if ( NUM(t) || (VR(t) != VR(p)) || ((vid)VR(p)->attr != V_PF) ) {
                    481:       *pr = t; return;
                    482:     } else {
                    483:       evalv(vl,VR(p),prec,&u); substr(vl,1,(Obj)t,VR(p),u,(Obj *)pr);
                    484:     }
                    485:   }
                    486: }
                    487:
                    488: void evalv(VL vl,V v,int prec,Obj *rp)
                    489: {
                    490:   PFINS ins,tins;
                    491:   PFAD ad,tad;
                    492:   PF pf;
                    493:   P t;
                    494:   int i;
                    495:
                    496:   if ( (vid)v->attr != V_PF ) {
                    497:     MKV(v,t); *rp = (Obj)t;
                    498:   } else {
                    499:     ins = (PFINS)v->priv; ad = ins->ad; pf = ins->pf;
                    500:     tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
                    501:     tins->pf = pf;
                    502:     for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) {
                    503:       tad[i].d = ad[i].d; evalr(vl,ad[i].arg,prec,&tad[i].arg);
                    504:     }
                    505:     evalins(tins,prec,rp);
                    506:   }
                    507: }
                    508:
                    509: void evalins(PFINS ins,int prec,Obj *rp)
                    510: {
                    511:   PF pf;
                    512:   PFINS tins;
                    513:   PFAD ad,tad;
                    514:   int i;
                    515:   Z q;
                    516:   V v;
                    517:   P x;
                    518:   NODE n0,n;
                    519:
                    520:   pf = ins->pf; ad = ins->ad;
                    521:   tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
                    522:   tins->pf = pf; tad = tins->ad;
                    523:   for ( i = 0; i < pf->argc; i++ ) {
                    524:     tad[i].d = ad[i].d; evalr(CO,ad[i].arg,prec,&tad[i].arg);
                    525:    }
                    526:   for ( i = 0; i < pf->argc; i++ )
                    527:     if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break;
                    528:   if ( (i != pf->argc) || !pf->pari ) {
                    529:     instoobj(tins,rp);
                    530:   } else {
                    531:     for ( n0 = 0, i = 0; i < pf->argc; i++ ) {
                    532:       NEXTNODE(n0,n); BDY(n) = (pointer)tad[i].arg;
                    533:     }
                    534:     if ( prec ) {
1.2       noro      535:       NEXTNODE(n0,n); STOZ(prec,q); BDY(n) = (pointer)q;
1.1       noro      536:     }
                    537:     if ( n0 )
                    538:       NEXT(n) = 0;
                    539:     (*pf->pari)(n0,rp);
                    540:   }
                    541: }
                    542:
                    543: void devalr(VL vl,Obj a,Obj *c)
                    544: {
                    545:   Obj nm,dn;
                    546:   double d;
                    547:   Real r,re,im;
                    548:   C z;
                    549:   int nid;
                    550:
                    551:   if ( !a )
                    552:     *c = 0;
                    553:   else {
                    554:     switch ( OID(a) ) {
                    555:       case O_N:
                    556:         nid = NID((Num)a);
                    557:         if ( nid == N_C ) {
                    558:           d = ToReal(((C)a)->r); MKReal(d,re);
                    559:           d = ToReal(((C)a)->i); MKReal(d,im);
                    560:           reimtocplx((Num)re,(Num)im,(Num *)&z);
                    561:           *c = (Obj)z;
                    562:         } else if ( nid == N_Q || nid == N_R || nid == N_B ) {
                    563:           d = ToReal(a);
                    564:           MKReal(d,r);
                    565:           *c = (Obj)r;
                    566:         } else
                    567:           error("devalr : unsupported");
                    568:         break;
                    569:       case O_P:
                    570:         devalp(vl,(P)a,(P *)c); break;
                    571:       case O_R:
                    572:         devalp(vl,NM((R)a),(P *)&nm);
                    573:         devalp(vl,DN((R)a),(P *)&dn);
                    574:         divr(vl,nm,dn,c);
                    575:         break;
                    576:       default:
                    577:         error("devalr : not implemented"); break;
                    578:     }
                    579:   }
                    580: }
                    581:
                    582: void devalp(VL vl,P p,P *pr)
                    583: {
                    584:   P t;
                    585:   DCP dc,dcr0,dcr;
                    586:   Obj u,s;
                    587:   double d;
                    588:   Real r;
                    589:
                    590:   if ( !p || NUM(p) ) {
                    591:     d = ToReal(p);
                    592:     MKReal(d,r);
                    593:     *pr = (P)r;
                    594:   } else {
                    595:     for ( dcr0 = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
                    596:       devalp(vl,COEF(dc),&t);
                    597:       if ( t ) {
                    598:         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
                    599:       }
                    600:     }
                    601:     if ( !dcr0 )
                    602:       *pr = 0;
                    603:     else {
                    604:       NEXT(dcr) = 0; MKP(VR(p),dcr0,t);
                    605:       if ( NUM(t) ) {
                    606:         d = ToReal((Num)t);
                    607:         MKReal(d,r);
                    608:         *pr = (P)r;
                    609:       } else if ( (VR(t) != VR(p)) || (VR(p)->attr != (pointer)V_PF) )
                    610:         *pr = t;
                    611:       else {
                    612:         devalv(vl,VR(p),&u);
                    613:         substr(vl,1,(Obj)t,VR(p),u,&s);
                    614:         if ( s && NUM(s) ) {
                    615:           d = ToReal((Num)s);
                    616:           MKReal(d,r);
                    617:           *pr = (P)r;
                    618:         } else
                    619:           *pr = (P)s;
                    620:       }
                    621:     }
                    622:   }
                    623: }
                    624:
                    625: void devalv(VL vl,V v,Obj *rp)
                    626: {
                    627:   PFINS ins,tins;
                    628:   PFAD ad,tad;
                    629:   PF pf;
                    630:   P t;
                    631:   int i;
                    632:
                    633:   if ( (vid)v->attr != V_PF ) {
                    634:     MKV(v,t); *rp = (Obj)t;
                    635:   } else {
                    636:     ins = (PFINS)v->priv; ad = ins->ad; pf = ins->pf;
                    637:     tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
                    638:     tins->pf = pf;
                    639:     for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) {
                    640:       tad[i].d = ad[i].d; devalr(vl,ad[i].arg,&tad[i].arg);
                    641:     }
                    642:     devalins(tins,rp);
                    643:   }
                    644: }
                    645:
                    646: void devalins(PFINS ins,Obj *rp)
                    647: {
                    648:   PFINS tins;
                    649:   PF pf;
                    650:   PFAD ad,tad;
                    651:   int i;
                    652:   Real r;
                    653:   double d;
                    654:   V v;
                    655:   P x;
                    656:
                    657:   pf = ins->pf; ad = ins->ad;
                    658:   tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
                    659:   tins->pf = pf; tad = tins->ad;
                    660:   for ( i = 0; i < pf->argc; i++ ) {
                    661:     tad[i].d = ad[i].d; devalr(CO,ad[i].arg,&tad[i].arg);
                    662:    }
                    663:   for ( i = 0; i < pf->argc; i++ )
                    664:     if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break;
                    665:   if ( (i != pf->argc) || !pf->libm ) {
                    666:     instoobj(tins,rp);
                    667:   } else {
                    668:     for ( i = 0; i < pf->argc; i++ )
                    669:       if ( tad[i].arg && NID((Num)tad[i].arg) == N_C )
                    670:         error("devalins : not supported");
                    671:     switch ( pf->argc ) {
                    672:       case 0:
                    673:         d = (*pf->libm)(); break;
                    674:       case 1:
                    675:         d = (*pf->libm)(ToReal(tad[0].arg)); break;
                    676:       case 2:
                    677:         d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg)); break;
                    678:       case 3:
                    679:         d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg),
                    680:           ToReal(tad[2].arg)); break;
                    681:       case 4:
                    682:         d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg),
                    683:           ToReal(tad[2].arg),ToReal(tad[3].arg)); break;
                    684:       default:
                    685:         error("devalins : not supported");
                    686:     }
                    687:     MKReal(d,r); *rp = (Obj)r;
                    688:   }
                    689: }
                    690:
                    691: extern int evalef,bigfloat;
                    692:
                    693: void simplify_elemfunc_ins(PFINS ins,Obj *rp)
                    694: {
                    695:   if ( evalef ) {
                    696:     if ( bigfloat ) evalins(ins,0,rp);
                    697:     else devalins(ins,rp);
                    698:   } else instoobj(ins,rp);
                    699: }
                    700:
                    701: void simplify_factorial_ins(PFINS ins,Obj *rp)
                    702: {
                    703:   PFAD ad;
                    704:   Obj a;
                    705:   Z q;
                    706:
                    707:   ad = ins->ad;
                    708:   a = ad[0].arg;
                    709:   if ( !ad[0].d && INT(a) && smallz((Z)a) ) {
1.2       noro      710:     factorialz(ZTOS((Z)a),&q);
1.1       noro      711:     *rp = (Obj)q;
                    712:   } else simplify_elemfunc_ins(ins,rp);
                    713: }
                    714:
                    715: void simplify_abs_ins(PFINS ins,Obj *rp)
                    716: {
                    717:   PFAD ad;
                    718:   Obj a;
                    719:   Q q;
                    720:   double t;
                    721:   Real r;
                    722:   struct oNODE arg0;
                    723:
                    724:   ad = ins->ad;
                    725:   a = ad[0].arg;
                    726:   if ( !ad[0].d && NUM(a) && (!a || RATN(a)) ) {
                    727:     if ( !a || sgnq((Q)a) > 0 ) *rp = (Obj)a;
                    728:     else {
                    729:       chsgnq((Q)a,&q); *rp = (Obj)q;
                    730:     }
                    731:   } else if ( !ad[0].d && REAL(a) ) {
                    732:     t = fabs(((Real)a)->body);
                    733:     MKReal(t,r); *rp = (Obj)r;
                    734:   } else if ( !ad[0].d && BIGFLOAT(a) ) {
                    735:     arg0.body = (pointer)a; arg0.next = 0;
1.4     ! noro      736:     mp_abs(&arg0,(Num *)rp);
1.3       kondoh    737: #if defined(INTERVAL)
                    738:   } else if ( !ad[0].d && ITVD(a) ) {
                    739:     absintvald((IntervalDouble)a,(IntervalDouble*)rp);
                    740:   } else if ( !ad[0].d && ITVF(a) ) {
                    741:     absintvalp((Itv)a,(Itv*)rp);
                    742: #endif
1.1       noro      743:   } else simplify_elemfunc_ins(ins,rp);
                    744: }
                    745:
                    746: void simplify_ins(PFINS ins,Obj *rp)
                    747: {
                    748:   V v;
                    749:   P t;
                    750:
                    751:   if ( ins->pf->simplify )
                    752:     (*ins->pf->simplify)(ins,rp);
                    753:   else {
                    754:     instoobj(ins,rp);
                    755:   }
                    756: }
                    757:
                    758: void instoobj(PFINS ins,Obj *rp)
                    759: {
                    760:   V v,newv;
                    761:   P t;
                    762:
                    763:   NEWV(v); NAME(v) = 0;
                    764:   v->attr = (pointer)V_PF; v->priv = (pointer)ins;
                    765:   appendpfins(v,&newv);
                    766:   MKV(newv,t);
                    767:   *rp = (Obj)t;
                    768: }
                    769:
                    770: void substfr(VL vl,Obj a,PF u,PF f,Obj *c)
                    771: {
                    772:   Obj nm,dn;
                    773:
                    774:   if ( !a )
                    775:     *c = 0;
                    776:   else {
                    777:     switch ( OID(a) ) {
                    778:       case O_N:
                    779:         *c = a; break;
                    780:       case O_P:
                    781:         substfp(vl,a,u,f,c); break;
                    782:       case O_R:
                    783:         substfp(vl,(Obj)NM((R)a),u,f,&nm); substfp(vl,(Obj)DN((R)a),u,f,&dn);
                    784:         divr(vl,nm,dn,c);
                    785:         break;
                    786:       default:
                    787:         error("substfr : not implemented"); break;
                    788:     }
                    789:   }
                    790: }
                    791:
                    792: void substfp(VL vl,Obj p,PF u,PF f,Obj *pr)
                    793: {
                    794:   V v;
                    795:   DCP dc;
                    796:   Obj a,c,m,s,t,p0;
                    797:   Z d;
                    798:   P x;
                    799:
                    800:   if ( !p )
                    801:     *pr = 0;
                    802:   else if ( NUM(p) )
                    803:     *pr = (Obj)p;
                    804:   else {
                    805:     v = VR((P)p); dc = DC((P)p);
                    806:     if ( (long)v->attr != V_PF ) {
                    807:       MKV(VR((P)p),x);
                    808:       for ( c = 0; dc; dc = NEXT(dc) ) {
                    809:         substfp(vl,(Obj)COEF(dc),u,f,&t);
                    810:         if ( DEG(dc) ) {
                    811:           gen_pwrr(vl,(Obj)x,(Obj)DEG(dc),&s);
                    812:           mulr(vl,s,t,&m);
                    813:           addr(vl,m,c,&a); c = a;
                    814:         } else {
                    815:           addr(vl,t,c,&a); c = a;
                    816:         }
                    817:       }
                    818:     } else {
                    819:       substfv(vl,v,u,f,&p0);
                    820:       substfp(vl,(Obj)COEF(dc),u,f,&c);
                    821:       for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc = NEXT(dc) ) {
                    822:           subz(d,DEG(dc),(Z *)&t);
                    823:           gen_pwrr(vl,p0,t,&s); mulr(vl,s,c,&m);
                    824:           substfp(vl,(Obj)COEF(dc),u,f,&t); addr(vl,m,t,&c);
                    825:       }
                    826:       if ( d ) {
                    827:         gen_pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m);
                    828:         c = m;
                    829:       }
                    830:     }
                    831:     *pr = c;
                    832:   }
                    833: }
                    834:
                    835: void substfv(VL vl,V v,PF u,PF f,Obj *c)
                    836: {
                    837:   P t;
                    838:   Obj r,s,w;
                    839:   int i,j;
                    840:   PFINS ins,tins;
                    841:   PFAD ad,tad;
                    842:
                    843:   ins = (PFINS)v->priv; ad = ins->ad;
                    844:   if ( ins->pf == u ) {
                    845:     if ( u->argc != f->argc )
                    846:       error("substfv : argument mismatch");
                    847:     if ( !f->body ) {
                    848:       mkpfins(f,f->args,&v); MKV(v,t); r = (Obj)t;
                    849:     } else
                    850:       r = f->body;
                    851:     for ( i = 0; i < f->argc; i++ )
                    852:       for ( j = 0; j < ad[i].d; j++ ) {
                    853:         derivr(vl,r,f->args[i],&s); r = s;
                    854:       }
                    855:     for ( i = 0; i < f->argc; i++ ) {
                    856:       substfr(vl,ad[i].arg,u,f,&w);
                    857:       substr(vl,0,r,f->args[i],w,&s); r = s;
                    858:     }
                    859:     *c = r;
                    860:   } else {
                    861:     tins = (PFINS)MALLOC(sizeof(PF)+f->argc*sizeof(struct oPFAD));
                    862:     tins->pf = ins->pf; tad = tins->ad;
                    863:     for ( i = 0; i < f->argc; i++ ) {
                    864:       tad[i].d = ad[i].d; substfr(vl,ad[i].arg,u,f,&tad[i].arg);
                    865:     }
                    866:     instoobj(tins,c);
                    867:   }
                    868: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>