[BACK]Return to f-itv.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Annotation of OpenXM_contrib2/asir2000/engine/f-itv.c, Revision 1.2

1.1       saito       1: /*
1.2     ! kondoh      2:  * $OpenXM: OpenXM_contrib2/asir2000/engine/f-itv.c,v 1.1 2000/12/22 10:03:28 saito Exp $
1.1       saito       3: */
                      4: #if defined(INTERVAL)
                      5: #include "ca.h"
                      6: #include "base.h"
                      7: #if PARI
                      8: #include "genpari.h"
                      9: #include "itv-pari.h"
                     10: extern long prec;
                     11: #endif
                     12:
                     13: void ToBf(Num a, BF *rp)
                     14: {
                     15:        GEN     pa, pb, pc;
                     16:        BF      bn, bd;
                     17:        BF      c;
                     18:        long    ltop, lbot;
                     19:
                     20:        if ( ! a ) {
                     21:                *rp = 0;
                     22:        }
                     23:        else switch ( NID(a) ) {
                     24:                case N_Q:
                     25:                        ltop = avma;
                     26:                        ritopa_i(NM((Q)a), SGN((Q)a), &pa);
                     27:                        pb = cgetr(prec);
                     28:                        mpaff(pa, pb);
                     29:                        if ( INT((Q)a) ) {
                     30:                                lbot = avma;
                     31:                                pb = gerepile(ltop, lbot, pb);
                     32:                                patori(pb, &c);
                     33:                                cgiv(pb);
                     34:                                *rp = c;
                     35:                        } else {
                     36:                                patori(pb, &bn);
                     37:                                ritopa_i(DN((Q)a), 1, &pa);
                     38:                                pb = cgetr(prec);
                     39:                                mpaff(pa, pb);
                     40:                                lbot = avma;
                     41:                                pb = gerepile(ltop, lbot, pb);
                     42:                                patori(pb, &bd);
                     43:                                cgiv(pb);
                     44:                                divbf((Num)bn,(Num)bd, (Num *)&c);
                     45:                                *rp = c;
                     46:                        }
                     47:                        break;
                     48:                case N_R:
                     49:                        pb = dbltor(BDY((Real)a));
                     50:                        patori(pb, &c);
                     51:                        cgiv(pb);
                     52:                        *rp = c;
                     53:                        break;
                     54:                case N_B:
                     55:                        *rp = (BF)a;
                     56:                        break;
                     57:                default:
                     58:                        error("ToBf : not implemented");
                     59:                        break;
                     60:        }
                     61: }
                     62:
                     63: void double2bf(double d, BF *rp)
                     64: {
                     65:        GEN     p;
                     66:
                     67:        p = dbltor(d);
                     68:        patori(p, rp);
                     69:        cgiv(p);
                     70: }
                     71:
                     72: double bf2double(BF a)
                     73: {
                     74:        GEN     p;
                     75:
                     76:        ritopa(a, &p);
                     77:        return rtodbl(p);
                     78: }
                     79:
                     80: void getulp(BF a, BF *au)
                     81: {
                     82:        GEN     b, c;
                     83:        long    lb, i;
                     84:
                     85:        ritopa(a, &b);
                     86:        lb = lg(b);
                     87:        c = cgetr(lb);
                     88:        c[1] = b[1];
                     89:        for (i=2;i<lb-1;i++) c[i] = 0;
                     90:        c[i] = 1;
                     91:        setsigne(c, 1);
                     92:        patori(c,au);
                     93:        cgiv(c);
                     94:        cgiv(b);
                     95: }
                     96:
1.2     ! kondoh     97: void addulp(IntervalBigFloat a, IntervalBigFloat *c)
1.1       saito      98: {
                     99:        Num     ai, as, aiu, asu, inf, sup;
                    100:
                    101:        itvtois((Itv)a,&ai,&as);
                    102:        if ( ai ) {
                    103:                getulp((BF)ai, (BF *)&aiu);
                    104:                subbf(ai,aiu,&inf);
                    105:        } else  inf = ai;
                    106:        if ( as ) {
                    107:                getulp((BF)as, (BF *)&asu);
                    108:                addbf(as,asu,&sup);
                    109:        } else  sup = as;
                    110:        istoitv(inf,sup, (Itv *)c);
                    111: }
                    112:
1.2     ! kondoh    113: void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1       saito     114: {
                    115:        Num     ai,as,bi,bs,mas,mbs,tmp;
                    116:        Num     inf,sup;
                    117:        GEN     pa, pb, z;
                    118:        long    ltop, lbot;
                    119:
                    120:        if ( !a )
                    121:                *c = b;
                    122:        else if ( !b )
                    123:                *c = a;
                    124:        else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
                    125:                addnum(0,(Num)a,(Num)b,(Num *)c);
                    126:        else {
                    127:                itvtois((Itv)a,&inf,&sup);
                    128:                ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
                    129:                itvtois((Itv)b,&inf,&sup);
                    130:                ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
                    131: #if 0
                    132: printexpr(CO, ai);
                    133: printexpr(CO, as);
                    134: printexpr(CO, bi);
                    135: printexpr(CO, bs);
                    136: #endif
                    137:
                    138: #if 1
                    139:                addnum(0,ai,bi,&inf);
                    140:                addnum(0,as,bs,&sup);
                    141:                istoitv(inf,sup,(Itv *)&tmp);
1.2     ! kondoh    142:                addulp((IntervalBigFloat)tmp, c);
1.1       saito     143:                return;
                    144: #else
                    145:                ltop = avma;
                    146:                ritopa(ai,&pa);
                    147:                ritopa(bi,&pb);
                    148:                lbot = avma;
                    149:                z = gerepile(ltop,lbot,PariAddDown(pa,pb));
                    150:                patori(z,&inf); cgiv(z);
                    151:
                    152:        /* MUST check if ai, as, bi, bs are bigfloat. */
                    153:
                    154:                /*  as + bs = ( - ( (-as) + (-bs) ) ) */
                    155:                chsgnbf(as,&mas);
                    156:                chsgnbf(bs,&mbs);
                    157:                ltop = avma;
                    158:                ritopa(mas,&pa);
                    159:                ritopa(mbs,&pb);
                    160:                lbot = avma;
                    161:                z = gerepile(ltop,lbot,PariAddDown(pa,pb));
                    162:                patori(z,&tmp); cgiv(z);
                    163:
                    164:                chsgnbf(tmp,&sup);
                    165:                istoitv(inf,sup,c);
                    166: #endif
                    167:        }
                    168: }
                    169:
1.2     ! kondoh    170: void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1       saito     171: {
                    172:        Num     ai,as,bi,bs,mas, mbs;
                    173:        Num     inf,sup,tmp;
                    174:        GEN     pa, pb, z;
                    175:        long    ltop, lbot;
                    176:
                    177:        if ( !a )
                    178:                chsgnnum((Num)b,(Num *)c);
                    179:        else if ( !b )
                    180:                *c = a;
                    181:        else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
                    182:                subnum(0,(Num)a,(Num)b,(Num *)c);
                    183:        else {
                    184:                itvtois((Itv)a,&inf,&sup);
                    185:                ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
                    186:                itvtois((Itv)b,&inf,&sup);
                    187:                ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
                    188: #if 1
                    189:                subnum(0,ai,bs,&inf);
                    190:                subnum(0,as,bi,&sup);
                    191:                istoitv(inf,sup,(Itv *)&tmp);
1.2     ! kondoh    192:                addulp((IntervalBigFloat)tmp, c);
1.1       saito     193: #else
                    194:
                    195: /* MUST check if ai, as, bi, bs are bigfloat. */
                    196:                /* ai - bs = ai + (-bs)  */
                    197:                chsgnbf(bs,&mbs);
                    198:                ltop = avma;
                    199:                ritopa(ai,&pa);
                    200:                ritopa(mbs,&pb);
                    201:                lbot = avma;
                    202:                z = gerepile(ltop,lbot,PariAddDown(pa,pb));
                    203:                patori(z,&inf); cgiv(z);
                    204:
                    205:                /* as - bi = ( - ( bi + (-as) ) ) */
                    206:                chsgnbf(as,&mas);
                    207:                ltop = avma;
                    208:                ritopa(mas,&pa);
                    209:                ritopa(bi,&pb);
                    210:                lbot = avma;
                    211:                z = gerepile(ltop,lbot,PariAddDown(pa,pb));
                    212:                patori(z,&tmp); cgiv(z);
                    213:
                    214:                chsgnbf(tmp,&sup);
                    215:                istoitv(inf,sup,c);
                    216: #endif
                    217:        }
                    218: }
                    219:
1.2     ! kondoh    220: void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1       saito     221: {
                    222:        Num     ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;
                    223:        Num     inf, sup;
                    224:        int     ba,bb;
                    225:
                    226:        if ( !a || !b )
                    227:                *c = 0;
                    228:        else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
                    229:                mulnum(0,(Num)a,(Num)b,(Num *)c);
                    230:        else {
                    231:                itvtois((Itv)a,&inf,&sup);
                    232:                ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
                    233:                itvtois((Itv)b,&inf,&sup);
                    234:                ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
                    235:
                    236:        /* MUST check if ai, as, bi, bs are bigfloat. */
                    237:                chsgnnum(as,&a2);
                    238:                chsgnnum(bs,&b2);
                    239:
                    240:
                    241:                if ( ba = (compnum(0,0,a2) > 0) ) {
                    242:                        a1 = ai;
                    243:                } else {
                    244:                        a1 = a2;
                    245:                        a2 = ai;
                    246:                }
                    247:                if ( bb = (compnum(0,0,b2) > 0) ) {
                    248:                        b1 = bi;
                    249:                } else {
                    250:                        b1 = b2;
                    251:                        b2 = bi;
                    252:                }
                    253:                mulnum(0,a2,b2,&t);
                    254:                subnum(0,0,t,&c2);
                    255:                if ( compnum(0,0,b1) > 0 ) {
                    256:                        mulnum(0,a2,b1,&t);
                    257:                        subnum(0,0,t,&c1);
                    258:                        if ( compnum(0,0,a1) > 0 ) {
                    259:                                mulnum(0,a1,b2,&t);
                    260:                                subnum(0,0,t,&p);
                    261:                                if ( compnum(0,c1,p) > 0 ) c1 = p;
                    262:                                mulnum(0,a1,b1,&t);
                    263:                                subnum(0,0,t,&p);
                    264:                                if ( compnum(0,c2,p) > 0 ) c2 = p;
                    265:                        }
                    266:                } else {
                    267:                        if ( compnum(0,0,a1) > 0 ) {
                    268:                                mulnum(0,a1,b2,&t);
                    269:                                subnum(0,0,t,&c1);
                    270:                        } else {
                    271:                                mulnum(0,a1,b1,&c1);
                    272:                        }
                    273:                }
                    274:                if ( ba == bb ) {
                    275:                        subnum(0,0,c2,&t);
                    276:                        istoitv(c1,t,(Itv *)&tmp);
                    277:                } else {
                    278:                        subnum(0,0,c1,&t);
                    279:                        istoitv(c2,t,(Itv *)&tmp);
                    280:                }
1.2     ! kondoh    281:                addulp((IntervalBigFloat)tmp, c);
1.1       saito     282:        }
                    283: }
                    284:
                    285: int     initvf(Num a, Itv b)
                    286: {
                    287:        Num inf, sup;
                    288:
                    289:        itvtois(b, &inf, &sup);
                    290:        if ( compnum(0,inf,a) > 0 ) return 0;
                    291:        else if ( compnum(0,a,sup) > 0 ) return 0;
                    292:        else return 1;
                    293: }
                    294:
                    295: int     itvinitvf(Itv a, Itv b)
                    296: {
                    297:        Num ai,as,bi,bs;
                    298:
                    299:        itvtois(a, &ai, &as);
                    300:        itvtois(b, &bi, &bs);
                    301:        if ( compnum(0,ai,bi) >= 0  && compnum(0,bs,as) >= 0 ) return 1;
                    302:        else return 0;
                    303: }
                    304:
1.2     ! kondoh    305: void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1       saito     306: {
                    307:        Num     ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;
                    308:        Num     inf, sup;
                    309:        int     ba,bb;
                    310:
                    311:        if ( !b )
                    312:                error("divitv : division by 0");
                    313:        else if ( !a )
                    314:                *c = 0;
                    315:        else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
                    316:                divnum(0,(Num)a,(Num)b,(Num *)c);
                    317:        else {
                    318:                itvtois((Itv)a,&inf,&sup);
                    319:                ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
                    320:                itvtois((Itv)b,&inf,&sup);
                    321:                ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
                    322: /* MUST check if ai, as, bi, bs are bigfloat. */
                    323:                chsgnnum(as,&a2);
                    324:                chsgnnum(bs,&b2);
                    325:                if ( ba = (compnum(0,0,a2) > 0) ) {
                    326:                        a1 = ai;
                    327:                } else {
                    328:                        a1 = a2;
                    329:                        a2 = ai;
                    330:                }
                    331:                if ( bb = (compnum(0,0,b2) > 0) ) {
                    332:                        b1 = bi;
                    333:                } else {
                    334:                        b1 = b2;
                    335:                        b2 = bi;
                    336:                }
                    337:                if ( compnum(0,0,b1) >= 0 ) {
                    338:                        error("divitv : division by interval including 0.");
                    339:                } else {
                    340:                        divnum(0,a2,b1,&c2);
                    341:                        if ( compnum(0,0,a1) > 0 ) {
                    342:                                divnum(0,a1,b1,&c1);
                    343:                        } else {
                    344:                                divnum(0,a1,b2,&t);
                    345:                                subnum(0,0,t,&c1);
                    346:                        }
                    347:                }
                    348:                if ( ba == bb ) {
                    349:                        subnum(0,0,c2,&t);
                    350:                        istoitv(c1,t,(Itv *)&tmp);
                    351:                } else {
                    352:                        subnum(0,0,c1,&t);
                    353:                        istoitv(c2,t,(Itv *)&tmp);
                    354:                }
1.2     ! kondoh    355:                addulp((IntervalBigFloat)tmp, c);
1.1       saito     356:        }
                    357: }
                    358:
                    359: void pwritvf(Itv a, Num e, Itv *c)
                    360: {
                    361:        int ei;
                    362:        Itv t;
                    363:
                    364:        if ( !e )
                    365:                *c = (Itv)ONE;
                    366:        else if ( !a )
                    367:                *c = 0;
                    368:        else if ( NID(a) <= N_B )
                    369:                pwrnum(0,(Num)a,e,(Num *)c);
                    370:        else if ( !INT(e) ) {
                    371: #if PARI && 0
                    372:                GEN pa,pe,z;
                    373:                int ltop,lbot;
                    374:
                    375:                ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
                    376:                z = gerepile(ltop,lbot,gpui(pa,pe,prec));
                    377:                patori(z,c); cgiv(z);
                    378: #else
                    379:                error("pwritv : can't calculate a fractional power");
                    380: #endif
                    381:        } else {
                    382:                ei = QTOS((Q)e);
                    383:                pwritv0f(a,ei,&t);
                    384:                if ( SGN((Q)e) < 0 )
                    385:                        divbf((Num)ONE,(Num)t,(Num *)c);
                    386:                else
                    387:                        *c = t;
                    388:        }
                    389: }
                    390:
                    391: void pwritv0f(Itv a, int e, Itv *c)
                    392: {
                    393:        Num inf, sup;
                    394:        Num ai,Xmin,Xmax;
1.2     ! kondoh    395:        IntervalBigFloat tmp;
1.1       saito     396:        Q       ne;
                    397:
                    398:        if ( e == 1 )
                    399:                *c = a;
                    400:        else {
                    401:                STOQ(e,ne);
                    402:                if ( !(e%2) ) {
                    403:                        if ( initvp(0,a) ) {
                    404:                                Xmin = 0;
                    405:                                chsgnnum(INF(a),&ai);
                    406:                                if ( compnum(0,ai,SUP(a)) > 0 ) {
                    407:                                        Xmax = ai;
                    408:                                } else {
                    409:                                        Xmax = SUP(a);
                    410:                                }
                    411:                        } else {
                    412:                                if ( compnum(0,INF(a),0) > 0 ) {
                    413:                                        Xmin = INF(a);
                    414:                                        Xmax = SUP(a);
                    415:                                } else {
                    416:                                        Xmin = SUP(a);
                    417:                                        Xmax = INF(a);
                    418:                                }
                    419:                        }
                    420:                } else {
                    421:                        Xmin = INF(a);
                    422:                        Xmax = SUP(a);
                    423:                }
                    424:                if ( ! Xmin )   inf = 0;
                    425:                else            pwrbf(Xmin,(Num)ne,&inf);
                    426:                if ( ! Xmax )   sup = 0;
                    427:                else            pwrbf(Xmax,(Num)ne,&sup);
                    428:                istoitv(inf,sup,(Itv *)&tmp);
1.2     ! kondoh    429:                addulp(tmp, (IntervalBigFloat *)c);
1.1       saito     430:        }
                    431: }
                    432:
                    433: void chsgnitvf(Itv a, Itv *c)
                    434: {
                    435:        Num inf,sup;
                    436:
                    437:        if ( !a )
                    438:                *c = 0;
                    439:        else if ( NID(a) <= N_B )
                    440:                chsgnnum((Num)a,(Num *)c);
                    441:        else {
                    442:                chsgnnum(INF((Itv)a),&inf);
                    443:                chsgnnum(SUP((Itv)a),&sup);
                    444:                istoitv(inf,sup,c);
                    445:        }
                    446: }
                    447:
                    448: int cmpitvf(Itv a, Itv b)
                    449: {
                    450:        Num     ai,as,bi,bs;
                    451:        int     s,t;
                    452:
                    453:        if ( !a ) {
                    454:                if ( !b || (NID(b)<=N_B) )
                    455:                        return compnum(0,(Num)a,(Num)b);
                    456:                else
                    457:                        return -1;
                    458:        } else if ( !b ) {
                    459:                if ( !a || (NID(a)<=N_B) )
                    460:                        return compnum(0,(Num)a,(Num)b);
                    461:                else
                    462:                        return initvp((Num)b,a);
                    463:        } else {
                    464:                itvtois(a,&ai,&as);
                    465:                itvtois(b,&bi,&bs);
1.2     ! kondoh    466:                s = compnum(0,ai,bs) ;
        !           467:                t = compnum(0,bi,as) ;
        !           468:                if ( s > 0 ) return 1;
        !           469:                else if ( t > 0 ) return -1;
        !           470:                else  return 0;
1.1       saito     471:        }
                    472: }
                    473:
                    474: void miditvf(Itv a, Num *b)
                    475: {
                    476:        Num     ai,as;
                    477:        Num     t;
                    478:        Q       TWO;
                    479:
                    480:        if ( ! a ) *b = 0;
                    481:        else if ( (NID(a) <= N_B) )
                    482:                *b = (Num)a;
                    483:        else {
                    484:                STOQ(2,TWO);
                    485:                itvtois(a,&ai,&as);
                    486:                addbf(ai,as,&t);
                    487:                divbf(t,(Num)TWO,b);
                    488:        }
                    489: }
                    490:
                    491: void cupitvf(Itv a, Itv b, Itv *c)
                    492: {
                    493:        Num     ai,as,bi,bs;
                    494:        Num     inf,sup;
                    495:
                    496:        itvtois(a,&ai,&as);
                    497:        itvtois(b,&bi,&bs);
                    498:        if ( compnum(0,ai,bi) > 0 )     inf = bi;
                    499:        else                            inf = ai;
                    500:        if ( compnum(0,as,bs) > 0 )     sup = as;
                    501:        else                            sup = bs;
                    502:        istoitv(inf,sup,c);
                    503: }
                    504:
                    505: void capitvf(Itv a, Itv b, Itv *c)
                    506: {
                    507:        Num     ai,as,bi,bs;
                    508:        Num     inf,sup;
                    509:
                    510:        itvtois(a,&ai,&as);
                    511:        itvtois(b,&bi,&bs);
                    512:        if ( compnum(0,ai,bi) > 0 )     inf = ai;
                    513:        else                            inf = bi;
                    514:        if ( compnum(0,as,bs) > 0 )     sup = bs;
                    515:        else                            sup = as;
                    516:        if ( compnum(0,inf,sup) > 0 )   *c = 0;
                    517:        else                            istoitv(inf,sup,c);
                    518: }
                    519:
                    520:
                    521: void widthitvf(Itv a, Num *b)
                    522: {
                    523:        Num     ai,as;
                    524:
                    525:        if ( ! a ) *b = 0;
                    526:        else if ( (NID(a) <= N_B) )
                    527:                *b = (Num)a;
                    528:        else {
                    529:                itvtois(a,&ai,&as);
                    530:                subnum(0,as,ai,b);
                    531:        }
                    532: }
                    533:
                    534: void absitvf(Itv a, Num *b)
                    535: {
                    536:        Num     ai,as,bi,bs;
                    537:
                    538:        if ( ! a ) *b = 0;
                    539:        else if ( (NID(a) <= N_B) ) {
                    540:                if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
                    541:                else *b = (Num)a;
                    542:        } else {
                    543:                itvtois(a,&ai,&as);
                    544:                if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
                    545:                else bi = ai;
                    546:                if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
                    547:                else bs = as;
                    548:                if ( compnum(0,bi,bs) > 0 )     *b = bi;
                    549:                else                            *b = bs;
                    550:        }
                    551: }
                    552:
                    553: void distanceitvf(Itv a, Itv b, Num *c)
                    554: {
                    555:        Num     ai,as,bi,bs;
                    556:        Num     di,ds;
                    557:        Itv     d;
                    558:
                    559:        itvtois(a,&ai,&as);
                    560:        itvtois(b,&bi,&bs);
                    561:        subnum(0,ai,bi,&di);
                    562:        subnum(0,as,bs,&ds);
                    563:        istoitv(di,ds,&d);
                    564:        absitvf(d,c);
                    565: }
                    566:
                    567: #endif

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