[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.9

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

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