[BACK]Return to isolv.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Annotation of OpenXM_contrib2/asir2000/builtin/isolv.c, Revision 1.4

1.1       saito       1: /*
1.4     ! saito       2:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/isolv.c,v 1.3 2003/10/23 01:32:59 saito Exp $
1.1       saito       3:  */
                      4:
                      5: #include "ca.h"
                      6: #include "parse.h"
                      7: #include "version.h"
                      8:
                      9: #if defined(INTERVAL)
                     10:
                     11: static void Solve(NODE, Obj *);
                     12: static void NSolve(NODE, Obj *);
                     13:
                     14: void Solve1(P, Q, pointer *);
                     15: void Sturm(P, VECT *);
                     16: void boundbody(P, Q *);
1.3       saito      17: void binary(int , MAT);
1.1       saito      18: void separate(Q, Q, VECT, Q, Q, int, int, MAT, int *);
                     19: void ueval(P, Q, Q *);
                     20:
                     21: struct ftab isolv_tab[] = {
                     22:        {"solve", Solve, 2},
                     23:        {"nsolve", NSolve, 2},
                     24:        {0,0,0},
                     25: };
                     26:
                     27: static void
                     28: Solve(arg, rp)
                     29: NODE arg;
                     30: Obj  *rp;
                     31: {
                     32:        pointer p, Eps;
                     33:        pointer    root;
                     34:        V       v;
                     35:        Q       eps;
                     36:
                     37:        p = (pointer)ARG0(arg);
                     38:        if ( !p ) {
                     39:                *rp = 0;
                     40:                return;
                     41:        }
                     42:        Eps = (pointer)ARG1(arg);
                     43:        asir_assert(Eps, O_N, "solve");
                     44:        if ( NID(Eps) != N_Q ) {
                     45:                fprintf(stderr,"solve arg 2 is required a rational number");
                     46:                error(" : invalid argument");
                     47:                return;
                     48:        }
                     49:        DUPQ((Q)Eps, eps);
                     50:        SGN(eps) = 1;
                     51:        switch (OID(p)) {
                     52:                case O_N:
                     53:                        *rp = 0;
                     54:                        break;
                     55:                case O_P:
                     56:                        Pvars(arg, &root);
                     57:                        if (NEXT(BDY((LIST)root)) != 0) {
                     58:                                fprintf(stderr,"solve arg 1 is univariate polynormial");
                     59:                                error(" : invalid argument");
                     60:                                break;
                     61:                        }
                     62:                        Solve1((P)p, eps, &root);
                     63:                        *rp = (Obj)root;
                     64:                        break;
                     65:                case O_LIST:
                     66:                        fprintf(stderr,"solve,");
                     67:                        error(" : Sorry, not yet implement of multivars");
                     68:                        break;
                     69:                defaults:
                     70:                        *rp = 0;
                     71:        }
                     72: }
                     73:
                     74: static void
                     75: NSolve(arg, rp)
                     76: NODE arg;
                     77: Obj  *rp;
                     78: {
1.3       saito      79:        pointer p, Eps;
1.4     ! saito      80:        pointer root;
        !            81:        LIST            listp;
1.3       saito      82:        V                       v;
                     83:        Q                       eps;
                     84:        NODE            n, n0, m0, m, ln0;
                     85:        Num             r;
                     86:        Itv             iv;
                     87:        BF                      breal;
1.1       saito      88:
                     89:        p = (pointer)ARG0(arg);
                     90:        if ( !p ) {
                     91:                *rp = 0;
                     92:                return;
                     93:        }
                     94:        Eps = (pointer)ARG1(arg);
                     95:        asir_assert(Eps, O_N, "solve");
                     96:        if ( NID(Eps) != N_Q ) {
                     97:                fprintf(stderr,"solve arg 2 is required a rational number");
                     98:                error(" : invalid argument");
                     99:                return;
                    100:        }
                    101:        DUPQ((Q)Eps, eps);
                    102:        SGN(eps) = 1;
                    103:        switch (OID(p)) {
                    104:                case O_N:
                    105:                        *rp = 0;
                    106:                        break;
                    107:                case O_P:
                    108:                        Pvars(arg, &root);
                    109:                        if (NEXT(BDY((LIST)root)) != 0) {
                    110:                                fprintf(stderr,"solve arg 1 is univariate polynormial");
                    111:                                error(" : invalid argument");
                    112:                                break;
                    113:                        }
                    114:                        Solve1((P)p, eps, &root);
                    115:                        for (m0 = BDY((LIST)root), n0 = 0; m0; m0 = NEXT(m0)) {
                    116:                                m = BDY((LIST)BDY(m0));
                    117:                                miditvp(BDY(m), &r);
                    118:                                ToBf(r, &breal);
                    119:                                NEXTNODE( n0, n );
                    120:                                MKNODE(ln0, breal, NEXT(m));
1.3       saito     121:                                MKLIST((LIST)listp, ln0);
1.1       saito     122:                                BDY(n) = (pointer)listp;
                    123:                        }
                    124:                        NEXT(n) = 0;
1.3       saito     125:                        MKLIST((LIST)listp,n0);
1.1       saito     126:                        *rp = (pointer)listp;
                    127:                        break;
                    128:                case O_LIST:
                    129:                        fprintf(stderr,"solve,");
                    130:                        error(" : Sorry, not yet implement of multivars");
                    131:                        break;
                    132:                defaults:
                    133:                        *rp = 0;
                    134:        }
                    135: }
                    136:
                    137: void
                    138: Solve1(inp, eps, rt)
                    139: P    inp;
                    140: Q    eps;
                    141: pointer *rt;
                    142: {
                    143:        P    p;
                    144:        Q    up, low, a;
                    145:        DCP fctp, onedeg, zerodeg;
                    146:        LIST listp;
                    147:        VECT sseq;
                    148:        MAT  root;
                    149:        int  chvu, chvl, pad, tnumb, numb, i, j;
                    150:        Itv  iv;
                    151:        NODE n0, n, ln0, ln;
                    152:
                    153:        boundbody(inp, &up);
                    154:        if (!up) {
                    155:                *rt = 0;
                    156:                return;
                    157:        }
                    158:        Sturm(inp, &sseq);
                    159:        DUPQ(up,low);
                    160:        SGN(low) = -1;
                    161:        chvu = stumq(sseq, up);
                    162:        chvl = stumq(sseq, low);
                    163:        tnumb = abs(chvu - chvl);
1.3       saito     164:        MKMAT(root, tnumb, 4);
1.1       saito     165:        pad = -1;
                    166:
                    167:        fctrp(CO,inp,&fctp);
                    168:        for (fctp = NEXT(fctp), i = 0; fctp; fctp = NEXT(fctp)) {
                    169:                p = COEF(fctp);
                    170:                onedeg = DC(p);
                    171:                if ( !cmpq(DEG(onedeg), ONE) ) {
                    172:                        pad++;
                    173:                        if ( !NEXT(onedeg) ) {
                    174:                                BDY(root)[pad][0] = 0;
                    175:                                BDY(root)[pad][1] = 0;
                    176:                                BDY(root)[pad][2] = DEG(fctp);
1.3       saito     177:                                BDY(root)[pad][3] = p;
1.1       saito     178:                        } else {
                    179:                                divq((Q)COEF(NEXT(onedeg)),(Q)COEF(onedeg),&a);
                    180:                                BDY(root)[pad][0] = a;
                    181:                                BDY(root)[pad][1] = BDY(root)[pad][0];
                    182:                                BDY(root)[pad][2] = DEG(fctp);
1.3       saito     183:                                BDY(root)[pad][3] = p;
1.1       saito     184:                        }
                    185:                        continue;
                    186:                }
                    187:                boundbody(p, &up);
                    188:                Sturm(p, &sseq);
                    189:                DUPQ(up,low);
                    190:                SGN(low) = -1;
                    191:                chvu = stumq(sseq, up);
                    192:                chvl = stumq(sseq, low);
                    193:                numb = abs(chvu - chvl);
                    194:                separate(DEG(fctp), eps, sseq, up, low, chvu, chvl, root, &pad);
                    195:        }
                    196:        for (i = 0; i < pad; i++) {
                    197:                for (j = i; j <= pad; j++) {
                    198:                        if (cmpq(BDY(root)[i][0], BDY(root)[j][0]) > 0) {
                    199:                                a = BDY(root)[i][0];
                    200:                                BDY(root)[i][0] = BDY(root)[j][0];
                    201:                                BDY(root)[j][0] = a;
                    202:                                a = BDY(root)[i][1];
                    203:                                BDY(root)[i][1] = BDY(root)[j][1];
                    204:                                BDY(root)[j][1] = a;
                    205:                                a = BDY(root)[i][2];
                    206:                                BDY(root)[i][2] = BDY(root)[j][2];
                    207:                                BDY(root)[j][2] = a;
1.3       saito     208:                                a = BDY(root)[i][3];
                    209:                                BDY(root)[i][3] = BDY(root)[j][3];
                    210:                                BDY(root)[j][3] = a;
                    211:                        }
                    212:                }
                    213:        }
                    214:        for (i = 0; i < pad; i++) {
                    215:                while(cmpq(BDY(root)[i][1], BDY(root)[i+1][0]) > 0 ) {
                    216:                        binary(i, root);
                    217:                        binary(i+1, root);
                    218:                        if ( cmpq(BDY(root)[i][0], BDY(root)[i+1][1]) > 0 ) {
                    219:                                a = BDY(root)[i][0];
                    220:                                BDY(root)[i][0] = BDY(root)[i+1][0];
                    221:                                BDY(root)[i+1][0] = a;
                    222:                                a = BDY(root)[i][1];
                    223:                                BDY(root)[i][1] = BDY(root)[i+1][1];
                    224:                                BDY(root)[i+1][1] = a;
                    225:                                a = BDY(root)[i][2];
                    226:                                BDY(root)[i][2] = BDY(root)[i+1][2];
                    227:                                BDY(root)[i+1][2] = a;
                    228:                                a = BDY(root)[i][3];
                    229:                                BDY(root)[i][3] = BDY(root)[i+1][3];
                    230:                                BDY(root)[i+1][3] = a;
                    231:                                break;
1.1       saito     232:                        }
                    233:                }
                    234:        }
                    235:        for (i = 0, n0 = 0; i <= pad; i++) {
                    236:                istoitv(BDY(root)[i][0], BDY(root)[i][1], &iv);
                    237:                NEXTNODE(n0,n);
                    238:                MKNODE(ln, BDY(root)[i][2], 0); MKNODE(ln0, iv, ln);
                    239:                MKLIST(listp, ln0);BDY(n) = (pointer)listp;
                    240:        }
                    241:        NEXT(n) = 0;
                    242:        MKLIST(listp,n0);
                    243:        *rt = (pointer)listp;
                    244: }
                    245:
                    246: void
                    247: separate(mult, eps, sseq, up, low, upn, lown, root, padp)
                    248: VECT sseq;
                    249: Q              mult, eps, up, low;
                    250: int    upn, lown;
                    251: MAT    root;
                    252: int    *padp;
                    253: {
                    254:        int de, midn;
1.3       saito     255:        Q   mid, e;
1.1       saito     256:        P   p;
                    257:
                    258:        de = abs(lown - upn);
                    259:        if (de == 0) return;
                    260:        if (de == 1) {
                    261:                (*padp)++;
                    262:                BDY(root)[*padp][0] = up;
                    263:                BDY(root)[*padp][1] = low;
1.3       saito     264:                BDY(root)[*padp][3] = (P *)sseq->body[0];
1.1       saito     265:                subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e );
                    266:                SGN(e) = 1;
                    267:                while (cmpq(e, eps) > 0) {
1.3       saito     268:                        binary(*padp, root);
1.1       saito     269:                        subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e);
                    270:                        SGN(e) = 1;
                    271:                }
                    272:                BDY(root)[*padp][2] = mult;
                    273:                return;
                    274:        }
                    275:        addq(up, low, &mid);
                    276:        divq(mid, TWO, &mid);
                    277:        midn = stumq(sseq, mid);
                    278:        separate(mult, eps, sseq, low, mid, lown, midn, root, padp);
                    279:        separate(mult, eps, sseq, mid, up, midn, upn, root, padp);
1.3       saito     280: }
                    281:
                    282: void
                    283: binary(indx, root)
                    284: int indx;
                    285: MAT root;
                    286: {
                    287:        Q       a, b, c, d, e;
                    288:        P       p;
1.4     ! saito     289:        p = (P)BDY(root)[indx][3];
1.3       saito     290:        addq(BDY(root)[indx][0], BDY(root)[indx][1], &c);
                    291:        divq(c, TWO, &d);
                    292:        ueval(p, BDY(root)[indx][1], &a);
                    293:        ueval(p, d, &b);
                    294:        if (SGN(a) == SGN(b)){
                    295:                BDY(root)[indx][1] = d;
                    296:        } else {
                    297:                BDY(root)[indx][0] = d;
                    298:        }
1.1       saito     299: }
                    300:
                    301: void
                    302: Sturm(p, ret)
                    303: P    p;
                    304: VECT *ret;
                    305: {
                    306:        P    g1,g2,q,r,s, *t;
                    307:        Q    a,b,c,d,h,l,m,x;
                    308:        V    v;
                    309:        VECT seq;
                    310:        int  i,j;
                    311:
                    312:        v = VR(p);
                    313:        t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P));
                    314:        g1 = t[0] = p; diffp(CO,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2;
                    315:        for ( i = 1, h = ONE, x = ONE; ; ) {
                    316:                if ( NUM(g2) ) break;
                    317:                subq(DEG(DC(g1)),DEG(DC(g2)),&d);
                    318:                l = (Q)LC(g2);
                    319:                if ( SGN(l) < 0 ) {
                    320:                        chsgnq(l,&a); l = a;
                    321:                }
                    322:                addq(d,ONE,&a); pwrq(l,a,&b); mulp(CO,(P)b,g1,(P *)&a);
                    323:                divsrp(CO,(P)a,g2,&q,&r);
                    324:                if ( !r ) break;
                    325:                chsgnp(r,&s); r = s; i++;
                    326:                if ( NUM(r) ) {
                    327:                        t[i] = r; break;
                    328:                }
                    329:                pwrq(h,d,&m); g1 = g2;
                    330:                mulq(m,x,&a); divsp(CO,r,(P)a,&g2); t[i] = g2;
                    331:                x = (Q)LC(g1);
                    332:                if ( SGN(x) < 0 ) {
                    333:                        chsgnq(x,&a); x = a;
                    334:                }
                    335:                pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);
                    336:        }
                    337:        MKVECT(seq,i+1);
                    338:        for ( j = 0; j <= i; j++ ) seq->body[j] = (pointer)t[j];
                    339:        *ret = seq;
                    340: }
                    341:
                    342: int
                    343: stumq(s, val)
                    344: VECT s;
                    345: Q val;
                    346: {
                    347:        int len, i, j, c;
                    348:        P   *ss;
                    349:        Q   a, a0;
                    350:
                    351:        len = s->len;
                    352:        ss = (P *)s->body;
                    353:        for ( j = 0; j < len; j++ ){
                    354:                ueval(ss[j],val,&a0);
                    355:                if (a0) break;
                    356:        }
                    357:        for ( i = j++, c =0; i < len; i++) {
                    358:                ueval( ss[i], val, &a);
                    359:                if ( a ) {
                    360:                        if( (SGN(a) > 0 && SGN(a0) < 0) || (SGN(a) < 0 && SGN(a0) > 0) ){
                    361:                                c++;
                    362:                                a0 = a;
                    363:                        }
                    364:                }
                    365:        }
                    366:        return c;
                    367: }
                    368:
                    369: void
                    370: boundbody(p, q)
                    371: P  p;
                    372: Q *q;
                    373: {
                    374:        Q               t, max, tmp;
                    375:        DCP     dc;
                    376:
                    377:        if ( !p )
                    378:                *q = 0;
                    379:        else if ( p->id == O_N )
                    380:                *q = 0;
                    381:        else {
                    382:                NEWQ(tmp);
                    383:                SGN(tmp)=1;
                    384:                for ( dc = DC(p), max=0; dc; dc = NEXT(dc) ) {
                    385:                        t = (Q)COEF(dc);
                    386:                        NM(tmp)=NM(t);
                    387:                        DN(tmp)=DN(t);
                    388:                        if ( cmpq(tmp, max) > 0 ) DUPQ(tmp, max);
                    389:                }
                    390:                addq(ONE, max, q);
                    391:        }
                    392: }
                    393:
                    394: void
                    395: ueval(p, q, rp)
                    396: P p;
                    397: Q q;
                    398: Q *rp;
                    399: {
                    400:        Q   d, d1, a, b, t;
                    401:        Q   deg, da;
                    402:        Q   nm, dn;
                    403:        DCP dc;
                    404:
                    405:        if ( !p ) *rp = 0;
                    406:        else if ( NUM(p) ) *rp = (Q)p;
                    407:        else {
                    408:                if ( q ) {
                    409:                        NTOQ( DN(q), 1, dn );
                    410:                        NTOQ( NM(q), SGN(q), nm );
                    411:                } else {
                    412:                        dn = 0;
                    413:                        nm = 0;
                    414:                }
                    415:                if ( !dn ) {
                    416:                        dc = DC(p); t = (Q)COEF(dc);
                    417:                        for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {
                    418:                                subq(d, DEG(dc), &d1); pwrq(nm, d1, &a);
                    419:                                mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);
                    420:                        }
                    421:                        if ( d ) {
                    422:                                pwrq(nm,d,&a); mulq(t,a,&b); t = b;
                    423:                        }
                    424:                        *rp = t;
                    425:                } else {
                    426:                        dc = DC(p); t = (Q)COEF(dc);
                    427:                        for ( d=deg= DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {
                    428:                                subq(d, DEG(dc), &d1); pwrq(nm, d1, &a);
                    429:                                mulq(t,a,&b);
                    430:                                subq(deg, DEG(dc), &d1); pwrq(dn, d1, &a);
                    431:                                mulq(a, (Q)COEF(dc), &da);
                    432:                                addq(b,da,&t);
                    433:                        }
                    434:                        if ( d ) {
                    435:                                pwrq(nm,d,&a); mulq(t,a,&b); t = b;
                    436:                        }
                    437:                        *rp = t;
                    438:                }
                    439:        }
                    440: }
                    441: #endif

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