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

Annotation of OpenXM_contrib2/asir2000/engine/d-itv.c, Revision 1.7

1.1       saito       1: /*
1.7     ! ohara       2:  * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.6 2015/08/14 13:51:54 fujimoto Exp $
1.1       saito       3: */
                      4: #if defined(INTERVAL)
                      5: #include <float.h>
                      6: #include "ca.h"
                      7: #include "base.h"
1.3       ohara       8: #if defined(PARI)
1.1       saito       9: #include "genpari.h"
                     10: #endif
                     11:
                     12: #if defined(ITVDEBUG)
                     13: void printbinint(int d)
                     14: {
                     15:        int     i, j, mask;
                     16:        union {
                     17:                int     x;
                     18:                char    c[sizeof(int)];
                     19:        } a;
                     20:
                     21:        a.x = d;
                     22:        for(i=sizeof(int)-1;i>=0;i--) {
                     23:                mask = 0x80;
                     24:                for(j=0;j<8;j++) {
                     25:                        if (a.c[i] & mask) fprintf(stderr,"1");
                     26:                        else fprintf(stderr,"0");
                     27:                        mask >>= 1;
                     28:                }
                     29:        }
                     30:        fprintf(stderr,"\n");
                     31: }
                     32: #endif
                     33:
                     34: double NatToRealUp(N a, int *expo)
                     35: {
                     36:        int *p;
                     37:        int l,all,i,j,s,tb,top,tail;
                     38:        unsigned int t,m[2];
                     39:        N       b,c;
                     40:        int     kk, carry, rem;
                     41:
                     42:        p = BD(a); l = PL(a) - 1;
                     43:        for ( top = 0, t = p[l]; t; t >>= 1, top++ );
                     44:        all = top + BSH*l; tail = (53-top)%BSH; j = l-(53-top)/BSH-1;
                     45:
                     46:
                     47:        if ( j >= 0 ) {
                     48:
                     49: #if defined(ITVDEBUG)
                     50:        fprintf(stderr," %d : tail = %d\n", j, tail);
                     51:        printbinint(p[j]);
                     52: #endif
                     53:                kk = (1<< (BSH - tail)) - 1;
                     54:                rem = p[j] & kk;
                     55:                if ( !rem )
                     56:                        for (i=0;i<j;i++)
                     57:                                if ( p[j] ) {
                     58:                                        carry = 1;
                     59:                                        break;
                     60:                                }
                     61:                else carry = 1;
                     62:                if ( carry ) {
                     63:                        b = NALLOC(j+1+1);
                     64:                        PL(b) = j+1;
                     65:                        p = BD(b);
                     66:                        for(i=0;i<j;i++) p[i] = 0;
                     67:                        p[j]=(1<< (BSH - tail));
                     68:
                     69:                        addn(a,b,&c);
                     70:
                     71:                        p = BD(c); l = PL(c) - 1;
                     72:                        for ( top = 0, t = p[l]; t; t >>= 1, top++ );
                     73:                        all = top + BSH*l; tail = (53-top)%BSH; i = l-(53-top)/BSH-1;
                     74:                } else i = j;
                     75:        } else i = j;
                     76:
                     77:
                     78:        m[1] = i < 0 ? 0 : p[i]>>(BSH-tail);
                     79:        for ( j = 1, i++, tb = tail; i <= l; i++ ) {
                     80:                s = 32-tb; t = i < 0 ? 0 : p[i];
                     81:                if ( BSH > s ) {
                     82:                        m[j] |= ((t&((1<<s)-1))<<tb);
                     83:                        if ( !j )
                     84:                                break;
                     85:                        else {
                     86:                                j--; m[j] = t>>s; tb = BSH-s;
                     87:                        }
                     88:                } else {
                     89:                        m[j] |= (t<<tb); tb += BSH;
                     90:                }
                     91:        }
                     92:        s = (all-1)+1023;
                     93:        m[0] = (m[0]&((1<<20)-1))|(MIN(2046,s)<<20); *expo = MAX(s-2046,0);
                     94: #ifdef vax
                     95:        t = m[0]; m[0] = m[1]; m[1] = t; itod(m);
                     96: #endif
1.7     ! ohara      97: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
1.1       saito      98:        t = m[0]; m[0] = m[1]; m[1] = t;
                     99: #endif
                    100:        return *((double *)m);
                    101: }
                    102:
                    103: static double  Q2doubleDown(Q a)
                    104: {
                    105:        double nm,dn,t,snm,rd;
                    106:        int enm,edn,e;
                    107:        unsigned int *p,s;
                    108:
                    109:        nm = (SGN(a)>0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm);
                    110:
                    111:        if ( INT(a) )
                    112:                if ( enm >= 1 )
                    113:                        error("Q2doubleDown : Overflow");
                    114:                else
                    115:                        return SGN(a)>0 ? nm : -nm;
                    116:        else {
                    117:                dn = (SGN(a)>0)? NatToReal(DN(a),&edn): NatToRealUp(DN(a),&edn);
                    118:
                    119:                if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
                    120:                        error("Q2doubleDown : Overflow"); /* XXX */
                    121:                else {
                    122:                        t = 0.0; p = (unsigned int *)&t;        /* Machine */
                    123:                        *p = ((e+1023)<<20);
                    124: #ifdef vax
                    125:                        s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
                    126: #endif
1.7     ! ohara     127: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) | defined(ANDROID)
1.1       saito     128:                        s = p[0]; p[0] = p[1]; p[1] = s;
                    129: #endif
                    130:                        FPMINUSINF
                    131:                        snm = (SGN(a)>0) ? nm : -nm;
                    132:                        rd  = snm / dn * t;
                    133:                        FPNEAREST
                    134:                        return rd;
                    135:                }
                    136:        }
                    137: }
                    138:
                    139:
                    140: static double  Q2doubleUp(Q a)
                    141: {
                    142:        double nm,dn,t,snm,rd;
                    143:        int enm,edn,e;
                    144:        unsigned int *p,s;
                    145:
                    146:        nm = (SGN(a)>0)? NatToRealUp(NM(a),&enm): NatToReal(NM(a),&enm);
                    147:
                    148:        if ( INT(a) )
                    149:                if ( enm >= 1 )
                    150:                        error("Q2doubleUp : Overflow");
                    151:                else
                    152:                        return SGN(a)>0 ? nm : -nm;
                    153:        else {
                    154:                dn = (SGN(a)>0)? NatToRealUp(DN(a),&edn): NatToReal(DN(a),&edn);
                    155:
                    156:                if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
                    157:                        error("Q2doubleUp : Overflow"); /* XXX */
                    158:                else {
                    159:                        t = 0.0; p = (unsigned int *)&t;        /* Machine */
                    160:                        *p = ((e+1023)<<20);
                    161: #ifdef vax
                    162:                        s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
                    163: #endif
1.7     ! ohara     164: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
1.1       saito     165:                        s = p[0]; p[0] = p[1]; p[1] = s;
                    166: #endif
                    167: #if 0
                    168:                        FPPLUSINF
                    169:                        snm = (SGN(a)>0) ? nm : -nm;
                    170:                        rd  = snm / dn * t;
                    171: #endif
                    172:
                    173:                        FPMINUSINF
                    174:                        snm = (SGN(a)>0) ? -nm : nm;
                    175:                        rd  = snm / dn * t;
                    176:                        FPNEAREST
                    177:                        return (-rd);
                    178:                }
                    179:        }
                    180: }
                    181:
                    182: static double  PARI2doubleDown(BF a)
                    183: {
                    184:         GEN     p;
                    185:        double  d;
                    186:
                    187:         ritopa(a, &p);
                    188:         d = rtodbl(p);
                    189:        cgiv(p);
                    190:        return d;
                    191: }
                    192:
                    193: static double  PARI2doubleUp(BF a)
                    194: {
                    195:        return PARI2doubleDown(a);
                    196: }
                    197:
                    198: double subulpd(double d)
                    199: {
                    200:        double inf;
                    201:
                    202:        FPMINUSINF
                    203:        inf = d - DBL_MIN;
                    204:        FPNEAREST
                    205:        return inf;
                    206: }
                    207:
                    208: double addulpd(double d)
                    209: {
                    210:        double md, sup;
                    211:
                    212:        md = -d;
                    213:        FPMINUSINF
                    214:        sup = md - DBL_MIN;
                    215:        FPNEAREST
                    216:        return (-sup);
                    217: }
                    218:
                    219: double ToRealDown(Num a)
                    220: {
                    221:        double inf;
                    222:
                    223:        if ( ! a ) return 0.0;
                    224:        switch ( NID(a) ) {
                    225:                case N_Q:
                    226:                        inf = Q2doubleDown((Q)a); break;
                    227:                case N_R:
                    228:                        inf = subulpd(BDY((Real)a)); break;
                    229:                case N_B:
                    230:                        inf = PARI2doubleDown((BF)a); break;
                    231:                case N_IP:
                    232:                        inf = ToRealDown(INF((Itv)a));
                    233:                        break;
1.2       kondoh    234:                case N_IntervalDouble:
                    235:                        inf = INF((IntervalDouble)a); break;
1.1       saito     236:                case N_A:
                    237:                default:
                    238:                        inf = 0.0;
                    239:                        error("ToRealDown: not supported operands.");
                    240:                        break;
                    241:        }
                    242:        return inf;
                    243: }
                    244:
                    245: double ToRealUp(Num a)
                    246: {
                    247:        double sup;
                    248:
                    249:        if ( ! a ) return 0.0;
                    250:        switch ( NID(a) ) {
                    251:                case N_Q:
                    252:                        sup = Q2doubleUp((Q)a); break;
                    253:                case N_R:
                    254:                        sup = addulpd(BDY((Real)a)); break;
                    255:                case N_B:
                    256:                        sup = PARI2doubleUp((BF)a); break;
                    257:                case N_IP:
                    258:                        sup = ToRealUp(SUP((Itv)a)); break;
1.2       kondoh    259:                case N_IntervalDouble:
                    260:                        sup = SUP((IntervalDouble)a); break;
1.1       saito     261:                case N_A:
                    262:                default:
                    263:                        sup = 0.0;
                    264:                        error("ToRealUp: not supported operands.");
                    265:                        break;
                    266:        }
                    267:        return sup;
                    268: }
                    269:
                    270:
                    271: void   Num2double(Num a, double *inf, double *sup)
                    272: {
                    273:        switch ( NID(a) ) {
                    274:                case N_Q:
                    275:                        *inf = Q2doubleDown((Q)a);
                    276:                        *sup = Q2doubleUp((Q)a);
                    277:                        break;
                    278:                case N_R:
                    279:                        *inf = BDY((Real)a);
                    280:                        *sup = BDY((Real)a);
                    281:                        break;
                    282:                case N_B:
                    283:                        *inf = PARI2doubleDown((BF)a);
                    284:                        *sup = PARI2doubleUp((BF)a);
                    285:                        break;
                    286:                case N_IP:
                    287:                        *inf = ToRealDown(INF((Itv)a));
                    288:                        *sup = ToRealUp(SUP((Itv)a));
                    289:                        break;
1.2       kondoh    290:                case N_IntervalDouble:
                    291:                        *inf = INF((IntervalDouble)a);
                    292:                        *sup = SUP((IntervalDouble)a);
1.1       saito     293:                        break;
                    294:                case N_A:
                    295:                default:
                    296:                        *inf = 0.0;
                    297:                        *sup = 0.0;
                    298:                        error("Num2double: not supported operands.");
                    299:                        break;
                    300:        }
                    301: }
                    302:
                    303:
1.2       kondoh    304: void additvd(Num a, Num b, IntervalDouble *c)
1.1       saito     305: {
                    306:        double  ai,as,mas, bi,bs;
                    307:        double  inf,sup;
                    308:
                    309:        if ( !a ) {
1.2       kondoh    310:                *c = (IntervalDouble)b;
1.1       saito     311:        } else if ( !b ) {
1.2       kondoh    312:                *c = (IntervalDouble)a;
1.1       saito     313: #if    0
                    314:        } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
                    315:                addnum(0,a,b,c);
                    316: #endif
                    317:        } else {
                    318:
                    319:                Num2double((Num)a,&ai,&as);
                    320:                Num2double((Num)b,&bi,&bs);
                    321:                mas = -as;
                    322:                FPMINUSINF
                    323:                inf = ai + bi;
                    324:                sup = mas - bs;         /*  as + bs = ( - ( (-as) - bs ) ) */
                    325:                FPNEAREST
1.2       kondoh    326:                MKIntervalDouble(inf,(-sup),*c);
1.1       saito     327:        }
                    328: }
                    329:
1.2       kondoh    330: void subitvd(Num a, Num b, IntervalDouble *c)
1.1       saito     331: {
                    332:        double  ai,as,mas, bi,bs;
                    333:        double  inf,sup;
                    334:
                    335:        if ( !a ) {
1.2       kondoh    336:                *c = (IntervalDouble)b;
1.1       saito     337:        } else if ( !b ) {
1.2       kondoh    338:                *c = (IntervalDouble)a;
1.1       saito     339: #if    0
                    340:        } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
                    341:                addnum(0,a,b,c);
                    342: #endif
                    343:        } else {
                    344:                Num2double(a,&ai,&as);
                    345:                Num2double(b,&bi,&bs);
                    346:                FPMINUSINF
                    347:                inf = ai - bs;
                    348:                sup = ( bi - as );      /* as - bi = ( - ( bi - as ) ) */
                    349:                FPNEAREST
1.2       kondoh    350:                MKIntervalDouble(inf,(-sup),*c);
1.1       saito     351:        }
                    352: }
                    353:
1.2       kondoh    354: void   mulitvd(Num a, Num b, IntervalDouble *c)
1.1       saito     355: {
                    356:        double  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p;
                    357:        double  inf, sup;
                    358:        int     ba,bb;
                    359:
                    360:        if ( !a || !b ) {
                    361:                *c = 0;
                    362: #if    0
                    363:        } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
                    364:                mulnum(0,a,b,c);
                    365: #endif
                    366:        } else {
                    367:                Num2double(a,&ai,&as);
                    368:                Num2double(b,&bi,&bs);
                    369:
                    370:                FPMINUSINF
                    371:
                    372:                a2 = -as;
                    373:                b2 = -bs;
                    374:
                    375:                if ( ba = ( a2 < 0.0 ) ) {
                    376:                        a1 = ai;
                    377:                } else {
                    378:                        a1 = a2;
                    379:                        a2 = ai;
                    380:                }
                    381:                if ( bb = ( b2 < 0.0 ) ) {
                    382:                        b1 = bi;
                    383:                } else {
                    384:                        b1 = b2;
                    385:                        b2 = bi;
                    386:                }
                    387:
                    388:                c2 = - a2 * b2;
                    389:                if ( b1 < 0.0 ) {
                    390:                        c1 = - a2 * b1;
                    391:                        if ( a1 < 0.0 ) {
                    392:                                p = - a1 * b2;
                    393:                                if ( p < c1 ) c1 = p;
                    394:                                p = - a1 * b1;
                    395:                                if ( p < c2 ) c2 = p;
                    396:                        }
                    397:                } else {
                    398:                        c1 = (a1<0.0)?( - a1 * b2 ):( a1 * b1 );
                    399:                }
                    400:                if ( ba == bb ) {
                    401:                        inf = c1;
                    402:                        sup = - c2;
                    403:                } else {
                    404:                        inf = c2;
                    405:                        sup = - c1;
                    406:                }
                    407:                FPNEAREST
1.2       kondoh    408:                MKIntervalDouble(inf,sup,*c);
1.1       saito     409:        }
                    410: }
                    411:
1.2       kondoh    412: int     initvd(Num a, IntervalDouble b)
1.1       saito     413: {
                    414:        Real inf, sup;
                    415:
1.2       kondoh    416:        if ( NID(b) == N_IntervalDouble ) {
1.1       saito     417:                MKReal(INF(b), inf);
                    418:                MKReal(SUP(b), sup);
                    419:        } else return 0;
                    420:        if ( compnum(0,(Num)inf,a) > 0 ) return 0;
                    421:        else if ( compnum(0,a,(Num)sup) > 0 ) return 0;
                    422:        else return 1;
                    423: }
                    424:
1.2       kondoh    425: void   divitvd(Num a, Num b, IntervalDouble *c)
1.1       saito     426: {
                    427:        double  ai,as,bi,bs,a1,a2,b1,b2,c1,c2;
                    428:        double  inf, sup;
                    429:        int     ba,bb;
                    430:
                    431:        if ( !b ) {
                    432:                *c = 0;
                    433:                error("divitvd : division by 0");
                    434:        } else if ( !a ) {
                    435:                *c = 0;
                    436: #if    0
                    437:        } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
                    438:                divnum(0,a,b,c);
                    439: #endif
                    440:        } else {
                    441:                Num2double(a,&ai,&as);
                    442:                Num2double(b,&bi,&bs);
                    443:
                    444:                FPMINUSINF
                    445:
                    446:                a2 = -as;
                    447:                b2 = -bs;
                    448:
                    449:                if ( ba = ( a2 < 0.0 ) ) {
                    450:                        a1 = ai;
                    451:                } else {
                    452:                        a1 = a2;
                    453:                        a2 = ai;
                    454:                }
                    455:                if ( bb = ( b2 < 0.0 ) ) {
                    456:                        b1 = bi;
                    457:                } else {
                    458:                        b1 = b2;
                    459:                        b2 = bi;
                    460:                }
                    461:
                    462:                if ( b1 <= 0.0 ) {
                    463:                        *c = 0;
                    464:                        error("divitvd : division by 0 interval");
                    465:                } else {
                    466:                        c2 = a2 / b1;
                    467:                        c1 = (a1<0.0)?( a1 / b1 ):( - a1 / b2 );
                    468:                }
                    469:                if ( ba == bb ) {
                    470:                        inf = c1;
                    471:                        sup = - c2;
                    472:                } else {
                    473:                        inf = c2;
                    474:                        sup = - c1;
                    475:                }
                    476:                FPNEAREST
1.2       kondoh    477:                MKIntervalDouble(inf,sup,*c);
1.1       saito     478:        }
                    479: }
                    480:
1.2       kondoh    481: void   pwritvd(Num a, Num e, IntervalDouble *c)
1.1       saito     482: {
                    483:        int     ei;
1.2       kondoh    484:        IntervalDouble  t;
1.1       saito     485:
                    486:        if ( !e ) {
1.2       kondoh    487:                *c = (IntervalDouble)ONE;
1.1       saito     488:        }  else if ( !a ) {
                    489:                *c = 0;
                    490: #if    0
                    491:        } else if ( NID(a) <= N_IP ) {
                    492:                pwrnum(0,a,e,c);
                    493: #endif
                    494:        } else if ( !INT(e) ) {
1.3       ohara     495: #if defined(PARI) && 0
1.4       ohara     496:                gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1       saito     497: #else
                    498:                error("pwritvd : can't calculate a fractional power");
                    499: #endif
                    500:        } else {
                    501:                ei = QTOS((Q)e);
1.2       kondoh    502:                pwritv0d((IntervalDouble)a,ei,&t);
1.1       saito     503:                if ( SGN((Q)e) < 0 )
                    504:                        divnum(0,(Num)ONE,(Num)t,(Num *)c);
                    505:                else
                    506:                        *c = t;
                    507:        }
                    508: }
                    509:
1.2       kondoh    510: void   pwritv0d(IntervalDouble a, int e, IntervalDouble *c)
1.1       saito     511: {
                    512:        double inf, sup;
                    513:        double t, Xmin, Xmax;
                    514:        int i;
                    515:
                    516:        if ( e == 1 )
                    517:                *c = a;
                    518:        else {
                    519:                if ( !(e%2) ) {
                    520:                        if ( initvd(0,a) ) {
                    521:                                Xmin = 0.0;
                    522:                                t = - INF(a);
                    523:                                Xmax = SUP(a);
                    524:                                if ( t > Xmax ) Xmax = t;
                    525:                        } else {
                    526:                                if ( INF(a) > 0.0 ) {
                    527:                                        Xmin = INF(a);
                    528:                                        Xmax = SUP(a);
                    529:                                } else {
                    530:                                        Xmin = - SUP(a);
                    531:                                        Xmax = - INF(a);
                    532:                                }
                    533:                        }
                    534:                } else {
                    535:                        Xmin = INF(a);
                    536:                        Xmax = SUP(a);
                    537:                }
                    538:                FPPLUSINF
                    539:                sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0;
                    540:                FPMINUSINF
                    541:                inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;
                    542:                FPNEAREST
1.2       kondoh    543:                MKIntervalDouble(inf,sup,*c);
1.1       saito     544:        }
                    545: }
                    546:
1.2       kondoh    547: void   chsgnitvd(IntervalDouble a, IntervalDouble *c)
1.1       saito     548: {
                    549:        double inf,sup;
                    550:
                    551:        if ( !a )
                    552:                *c = 0;
                    553: #if    0
                    554:        else if ( NID(a) <= N_IP )
                    555:                chsgnnum(a,c);
                    556: #endif
                    557:        else {
                    558:                inf = - SUP(a);
                    559:                sup = - INF(a);
1.2       kondoh    560:                MKIntervalDouble(inf,sup,*c);
1.1       saito     561:        }
                    562: }
                    563:
1.2       kondoh    564: int    cmpitvd(IntervalDouble a, IntervalDouble b)
                    565: /*    a > b:  1       */
                    566: /*    a = b:  0       */
                    567: /*    a < b: -1       */
1.1       saito     568: {
                    569:        double  ai,as,bi,bs;
                    570:
                    571:        if ( !a ) {
                    572:                if ( !b || (NID(b)<=N_IP) ) {
                    573:                        return compnum(0,(Num)a,(Num)b);
                    574:                } else {
                    575:                        bi = INF(b);
                    576:                        bs = SUP(b);
                    577:                        if ( bi > 0.0 ) return -1;
                    578:                        else if ( bs < 0.0 ) return 1;
                    579:                        else return 0;
                    580:                }
                    581:        } else if ( !b ) {
                    582:                if ( !a || (NID(a)<=N_IP) ) {
                    583:                        return compnum(0,(Num)a,(Num)b);
                    584:                } else {
                    585:                        ai = INF(a);
                    586:                        as = SUP(a);
                    587:                        if ( ai > 0.0 ) return 1;
                    588:                        else if ( as < 0.0 ) return -1;
                    589:                        else return 0;
                    590:                }
                    591:        } else {
                    592:                bi = INF(b);
                    593:                bs = SUP(b);
                    594:                ai = INF(a);
                    595:                as = SUP(a);
                    596:                if ( ai > bs ) return 1;
                    597:                else if ( bi > as ) return -1;
                    598:                else return 0;
                    599:        }
                    600: }
                    601:
1.2       kondoh    602: void   miditvd(IntervalDouble a, Num *b)
1.1       saito     603: {
                    604:        double  t;
                    605:        Real    rp;
                    606:
                    607:        if ( ! a ) *b = 0;
                    608:        else if ( (NID(a) <= N_IP) )
                    609:                *b = (Num)a;
                    610:        else {
                    611:                t = (INF(a) + SUP(a))/2.0;
                    612:                MKReal(t,rp);
                    613:                *b = (Num)rp;
                    614:        }
                    615: }
                    616:
1.2       kondoh    617: void   cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
1.1       saito     618: {
                    619:        double  ai,as,bi,bs;
                    620:        double  inf,sup;
                    621:
                    622:        ai = INF(a);
                    623:        as = SUP(a);
                    624:        bi = INF(b);
                    625:        bs = SUP(b);
                    626:        inf = MIN(ai,bi);
                    627:        sup = MAX(as,bs);
1.2       kondoh    628:        MKIntervalDouble(inf,sup,*c);
1.1       saito     629: }
                    630:
1.2       kondoh    631: void   capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
1.1       saito     632: {
                    633:        double  ai,as,bi,bs;
                    634:        double  inf,sup;
                    635:
                    636:        ai = INF(a);
                    637:        as = SUP(a);
                    638:        bi = INF(b);
                    639:        bs = SUP(b);
                    640:        inf = MAX(ai,bi);
                    641:        sup = MIN(as,bs);
                    642:        if ( inf > sup ) *c = 0;
1.2       kondoh    643:        else MKIntervalDouble(inf,sup,*c);
1.1       saito     644: }
                    645:
                    646:
1.2       kondoh    647: void   widthitvd(IntervalDouble a, Num *b)
1.1       saito     648: {
                    649:        double  t;
                    650:        Real    rp;
                    651:
                    652:        if ( ! a ) *b = 0;
                    653:        else {
                    654:                t = SUP(a) - INF(a);
                    655:                MKReal(t,rp);
                    656:                *b = (Num)rp;
                    657:        }
                    658: }
                    659:
1.2       kondoh    660: void   absitvd(IntervalDouble a, Num *b)
1.1       saito     661: {
                    662:        double  ai,as,bi,bs;
                    663:        double  t;
                    664:        Real    rp;
                    665:
                    666:        if ( ! a ) *b = 0;
                    667:        else {
                    668:                ai = INF(a);
                    669:                as = SUP(a);
                    670:                if (ai < 0) bi = -ai;   else    bi = ai;
                    671:                if (as < 0) bs = -as;   else    bs = as;
                    672:                if ( bi > bs )  t = bi; else    t = bs;
                    673:                MKReal(t,rp);
                    674:                *b = (Num)rp;
                    675:        }
                    676: }
                    677:
1.2       kondoh    678: void   distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)
1.1       saito     679: {
                    680:        double  ai,as,bi,bs;
                    681:        double  di,ds;
                    682:        double  t;
                    683:        Real    rp;
                    684:
                    685:        ai = INF(a);
                    686:        as = SUP(a);
                    687:        bi = INF(b);
                    688:        bs = SUP(b);
                    689:        di = bi - ai;
                    690:        ds = bs - as;
                    691:
                    692:        if (di < 0) di = -di;
                    693:        if (ds < 0) ds = -ds;
                    694:        if (di > ds)    t = di; else    t = ds;
                    695:        MKReal(t,rp);
                    696:        *c = (Num)rp;
                    697: }
                    698:
                    699: #endif

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