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

Annotation of OpenXM_contrib2/asir2018/engine/d-itv.c, Revision 1.1

1.1     ! noro        1: /*
        !             2:  * $OpenXM$
        !             3: */
        !             4: #if defined(INTERVAL)
        !             5: #include <float.h>
        !             6: #include "ca.h"
        !             7: #include "base.h"
        !             8: #if defined(PARI)
        !             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
        !            97: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
        !            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
        !           127: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) | defined(ANDROID)
        !           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
        !           164: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
        !           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;
        !           234:     case N_IntervalDouble:
        !           235:       inf = INF((IntervalDouble)a); break;
        !           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;
        !           259:     case N_IntervalDouble:
        !           260:       sup = SUP((IntervalDouble)a); break;
        !           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;
        !           290:     case N_IntervalDouble:
        !           291:       *inf = INF((IntervalDouble)a);
        !           292:       *sup = SUP((IntervalDouble)a);
        !           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:
        !           304: void additvd(Num a, Num b, IntervalDouble *c)
        !           305: {
        !           306:   double  ai,as,mas, bi,bs;
        !           307:   double  inf,sup;
        !           308:
        !           309:   if ( !a ) {
        !           310:     *c = (IntervalDouble)b;
        !           311:   } else if ( !b ) {
        !           312:     *c = (IntervalDouble)a;
        !           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
        !           326:     MKIntervalDouble(inf,(-sup),*c);
        !           327:   }
        !           328: }
        !           329:
        !           330: void subitvd(Num a, Num b, IntervalDouble *c)
        !           331: {
        !           332:   double  ai,as,mas, bi,bs;
        !           333:   double  inf,sup;
        !           334:
        !           335:   if ( !a ) {
        !           336:     *c = (IntervalDouble)b;
        !           337:   } else if ( !b ) {
        !           338:     *c = (IntervalDouble)a;
        !           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
        !           350:     MKIntervalDouble(inf,(-sup),*c);
        !           351:   }
        !           352: }
        !           353:
        !           354: void  mulitvd(Num a, Num b, IntervalDouble *c)
        !           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
        !           408:     MKIntervalDouble(inf,sup,*c);
        !           409:   }
        !           410: }
        !           411:
        !           412: int     initvd(Num a, IntervalDouble b)
        !           413: {
        !           414:   Real inf, sup;
        !           415:
        !           416:   if ( NID(b) == N_IntervalDouble ) {
        !           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:
        !           425: void  divitvd(Num a, Num b, IntervalDouble *c)
        !           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
        !           477:     MKIntervalDouble(inf,sup,*c);
        !           478:   }
        !           479: }
        !           480:
        !           481: void  pwritvd(Num a, Num e, IntervalDouble *c)
        !           482: {
        !           483:   int  ei;
        !           484:   IntervalDouble  t;
        !           485:
        !           486:   if ( !e ) {
        !           487:     *c = (IntervalDouble)ONE;
        !           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) ) {
        !           495: #if defined(PARI) && 0
        !           496:     gpui_ri((Obj)a,(Obj)c,(Obj *)c);
        !           497: #else
        !           498:     error("pwritvd : can't calculate a fractional power");
        !           499: #endif
        !           500:   } else {
        !           501:     ei = QTOS((Q)e);
        !           502:     pwritv0d((IntervalDouble)a,ei,&t);
        !           503:     if ( SGN((Q)e) < 0 )
        !           504:       divnum(0,(Num)ONE,(Num)t,(Num *)c);
        !           505:     else
        !           506:       *c = t;
        !           507:   }
        !           508: }
        !           509:
        !           510: void  pwritv0d(IntervalDouble a, int e, IntervalDouble *c)
        !           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
        !           543:     MKIntervalDouble(inf,sup,*c);
        !           544:   }
        !           545: }
        !           546:
        !           547: void  chsgnitvd(IntervalDouble a, IntervalDouble *c)
        !           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);
        !           560:     MKIntervalDouble(inf,sup,*c);
        !           561:   }
        !           562: }
        !           563:
        !           564: int  cmpitvd(IntervalDouble a, IntervalDouble b)
        !           565: /*    a > b:  1       */
        !           566: /*    a = b:  0       */
        !           567: /*    a < b: -1       */
        !           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:
        !           602: void  miditvd(IntervalDouble a, Num *b)
        !           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:
        !           617: void  cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
        !           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);
        !           628:   MKIntervalDouble(inf,sup,*c);
        !           629: }
        !           630:
        !           631: void  capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
        !           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;
        !           643:   else MKIntervalDouble(inf,sup,*c);
        !           644: }
        !           645:
        !           646:
        !           647: void  widthitvd(IntervalDouble a, Num *b)
        !           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:
        !           660: void  absitvd(IntervalDouble a, Num *b)
        !           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:
        !           678: void  distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)
        !           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>