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

Annotation of OpenXM_contrib2/asir2000/engine/t-itv.c, Revision 1.6

1.1       saito       1: /*
1.6     ! noro        2:  * $OpenXM: OpenXM_contrib2/asir2000/engine/t-itv.c,v 1.5 2009/03/27 14:42:29 ohara Exp $
1.1       saito       3: */
                      4: #if defined(INTERVAL)
                      5: #include "ca.h"
                      6: #include "base.h"
1.3       ohara       7: #if defined(PARI)
1.1       saito       8: #include "genpari.h"
                      9: #endif
                     10:
                     11: void itvtois(Itv a, Num *inf, Num *sup)
                     12: {
1.6     ! noro       13:   if ( !a )
        !            14:     *inf = *sup = (Num)0;
        !            15:   else if ( NID(a) <= N_B ) {
        !            16:     *inf = (Num)a; *sup = (Num)a;
        !            17:   } else {
        !            18:     *inf = INF(a);
        !            19:     *sup = SUP(a);
        !            20:   }
1.1       saito      21: }
                     22:
                     23: void istoitv(Num inf, Num sup, Itv *rp)
                     24: {
1.6     ! noro       25:   Itv c;
1.1       saito      26:
1.6     ! noro       27:   if ( !inf && !sup )
        !            28:     *rp = 0;
        !            29:   else {
        !            30:     NEWItvP(c);
        !            31:     if ( compnum(0,inf,sup) >= 0 ) {
        !            32:       INF(c) = sup; SUP(c) = inf;
        !            33:     } else {
        !            34:       INF(c) = inf; SUP(c) = sup;
        !            35:     }
        !            36:     *rp = c;
        !            37:   }
1.1       saito      38: }
                     39:
                     40: void additvp(Itv a, Itv b, Itv *c)
                     41: {
1.6     ! noro       42:   Num  ai,as,bi,bs;
        !            43:   Num  inf,sup;
1.1       saito      44:
1.6     ! noro       45:   if ( !a )
        !            46:     *c = b;
        !            47:   else if ( !b )
        !            48:     *c = a;
        !            49:   else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
        !            50:     addnum(0,a,b,c);
        !            51:   else if ( (NID(a) == N_IP) && (NID(b) == N_R )
        !            52:     ||(NID(a) == N_R ) && (NID(b) == N_IP) )
        !            53:     additvd((Num)a,(Num)b,(IntervalDouble *)c);
        !            54:   else {
        !            55:     itvtois(a,&ai,&as);
        !            56:     itvtois(b,&bi,&bs);
        !            57:     addnum(0,ai,bi,&inf);
        !            58:     addnum(0,as,bs,&sup);
        !            59:     istoitv(inf,sup,c);
        !            60:   }
1.1       saito      61: }
                     62:
                     63: void subitvp(Itv a, Itv b, Itv *c)
                     64: {
1.6     ! noro       65:   Num  ai,as,bi,bs;
        !            66:   Num  inf,sup;
1.1       saito      67:
1.6     ! noro       68:   if ( !a )
        !            69:     chsgnnum(b,c);
        !            70:   else if ( !b )
        !            71:     *c = a;
        !            72:   else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
        !            73:     subnum(0,a,b,c);
        !            74:   else if ( (NID(a) == N_IP) && (NID(b) == N_R )
        !            75:     ||(NID(a) == N_R ) && (NID(b) == N_IP) )
        !            76:     subitvd((Num)a,(Num)b,(IntervalDouble *)c);
        !            77:   else {
        !            78:     itvtois(a,&ai,&as);
        !            79:     itvtois(b,&bi,&bs);
        !            80:     subnum(0,ai,bs,&inf);
        !            81:     subnum(0,as,bi,&sup);
        !            82:     istoitv(inf,sup,c);
        !            83:   }
1.1       saito      84: }
                     85:
                     86: void mulitvp(Itv a, Itv b, Itv *c)
                     87: {
1.6     ! noro       88:   Num  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t;
        !            89:   int  ba,bb;
1.1       saito      90:
1.6     ! noro       91:   if ( !a || !b )
        !            92:     *c = 0;
        !            93:   else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
        !            94:     mulnum(0,a,b,c);
        !            95:   else if ( (NID(a) == N_IP) && (NID(b) == N_R )
        !            96:     ||(NID(a) == N_R ) && (NID(b) == N_IP) )
        !            97:     mulitvd((Num)a,(Num)b,(IntervalDouble *)c);
        !            98:   else {
        !            99:     itvtois(a,&ai,&as);
        !           100:     itvtois(b,&bi,&bs);
        !           101:     chsgnnum(as,&a2);
        !           102:     chsgnnum(bs,&b2);
        !           103:     if ( ba = (compnum(0,0,a2) > 0) ) {
        !           104:       a1 = ai;
        !           105:     } else {
        !           106:       a1 = a2;
        !           107:       a2 = ai;
        !           108:     }
        !           109:     if ( bb = (compnum(0,0,b2) > 0) ) {
        !           110:       b1 = bi;
        !           111:     } else {
        !           112:       b1 = b2;
        !           113:       b2 = bi;
        !           114:     }
        !           115:     mulnum(0,a2,b2,&t);
        !           116:     subnum(0,0,t,&c2);
        !           117:     if ( compnum(0,0,b1) > 0 ) {
        !           118:       mulnum(0,a2,b1,&t);
        !           119:       subnum(0,0,t,&c1);
        !           120:       if ( compnum(0,0,a1) > 0 ) {
        !           121:         mulnum(0,a1,b2,&t);
        !           122:         subnum(0,0,t,&p);
        !           123:         if ( compnum(0,c1,p) > 0 ) c1 = p;
        !           124:         mulnum(0,a1,b1,&t);
        !           125:         subnum(0,0,t,&p);
        !           126:         if ( compnum(0,c2,p) > 0 ) c2 = p;
        !           127:       }
        !           128:     } else {
        !           129:       if ( compnum(0,0,a1) > 0 ) {
        !           130:         mulnum(0,a1,b2,&t);
        !           131:         subnum(0,0,t,&c1);
        !           132:       } else {
        !           133:         mulnum(0,a1,b1,&c1);
        !           134:       }
        !           135:     }
        !           136:     if ( ba == bb ) {
        !           137:       subnum(0,0,c2,&t);
        !           138:       istoitv(c1,t,c);
        !           139:     } else {
        !           140:       subnum(0,0,c1,&t);
        !           141:       istoitv(c2,t,c);
        !           142:     }
        !           143:   }
1.1       saito     144: }
                    145:
                    146: int     initvp(Num a, Itv b)
                    147: {
1.6     ! noro      148:   Num inf, sup;
1.1       saito     149:
1.6     ! noro      150:   itvtois(b, &inf, &sup);
        !           151:   if ( compnum(0,inf,a) > 0 ) return 0;
        !           152:   else if ( compnum(0,a,sup) > 0 ) return 0;
        !           153:   else return 1;
1.1       saito     154: }
                    155:
                    156: int     itvinitvp(Itv a, Itv b)
                    157: {
1.6     ! noro      158:   Num ai,as,bi,bs;
1.1       saito     159:
1.6     ! noro      160:   itvtois(a, &ai, &as);
        !           161:   itvtois(b, &bi, &bs);
        !           162:   if ( compnum(0,ai,bi) >= 0  && compnum(0,bs,as) >= 0 ) return 1;
        !           163:   else return 0;
1.1       saito     164: }
                    165:
                    166: void divitvp(Itv a, Itv b, Itv *c)
                    167: {
1.6     ! noro      168:   Num  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t;
        !           169:   int  ba,bb;
1.1       saito     170:
1.6     ! noro      171:   if ( !b )
        !           172:     error("divitv : division by 0");
        !           173:   else if ( !a )
        !           174:     *c = 0;
        !           175:   else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
        !           176:     divnum(0,a,b,c);
        !           177:   else if ( (NID(a) == N_IP) && (NID(b) == N_R )
        !           178:     ||(NID(a) == N_R ) && (NID(b) == N_IP) )
        !           179:     divitvd((Num)a,(Num)b,(IntervalDouble *)c);
        !           180:   else {
        !           181:     itvtois(a,&ai,&as);
        !           182:     itvtois(b,&bi,&bs);
        !           183:     chsgnnum(as,&a2);
        !           184:     chsgnnum(bs,&b2);
        !           185:     if ( ba = (compnum(0,0,a2) > 0) ) {
        !           186:       a1 = ai;
        !           187:     } else {
        !           188:       a1 = a2;
        !           189:       a2 = ai;
        !           190:     }
        !           191:     if ( bb = (compnum(0,0,b2) > 0) ) {
        !           192:       b1 = bi;
        !           193:     } else {
        !           194:       b1 = b2;
        !           195:       b2 = bi;
        !           196:     }
        !           197:     if ( compnum(0,0,b1) >= 0 ) {
        !           198:       error("divitv : division by interval including 0.");
        !           199:     } else {
        !           200:       divnum(0,a2,b1,&c2);
        !           201:       if ( compnum(0,0,a1) > 0 ) {
        !           202:         divnum(0,a1,b1,&c1);
        !           203:       } else {
        !           204:         divnum(0,a1,b2,&t);
        !           205:         subnum(0,0,t,&c1);
        !           206:       }
        !           207:     }
        !           208:     if ( ba == bb ) {
        !           209:       subnum(0,0,c2,&t);
        !           210:       istoitv(c1,t,c);
        !           211:     } else {
        !           212:       subnum(0,0,c1,&t);
        !           213:       istoitv(c2,t,c);
        !           214:     }
        !           215:   }
1.1       saito     216: }
                    217:
                    218: void pwritvp(Itv a, Num e, Itv *c)
                    219: {
1.6     ! noro      220:   int ei;
        !           221:   Itv t;
1.1       saito     222:
1.6     ! noro      223:   if ( !e )
        !           224:     *c = (Itv)ONE;
        !           225:   else if ( !a )
        !           226:     *c = 0;
        !           227:   else if ( NID(a) <= N_B )
        !           228:     pwrnum(0,a,e,c);
        !           229:   else if ( !INT(e) ) {
1.3       ohara     230: #if defined(PARI) && 0
1.6     ! noro      231:     gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1       saito     232: #else
1.6     ! noro      233:     error("pwritv : can't calculate a fractional power");
1.1       saito     234: #endif
1.6     ! noro      235:   } else {
        !           236:     ei = QTOS((Q)e);
        !           237:     pwritv0p(a,ei,&t);
        !           238:     if ( SGN((Q)e) < 0 )
        !           239:       divnum(0,(Num)ONE,(Num)t,c);
        !           240:     else
        !           241:       *c = t;
        !           242:   }
1.1       saito     243: }
                    244:
                    245: void pwritv0p(Itv a, int e, Itv *c)
                    246: {
1.6     ! noro      247:   Num inf, sup;
        !           248:   Num ai,Xmin,Xmax;
        !           249:   Q  ne;
        !           250:
        !           251:   if ( e == 1 )
        !           252:     *c = a;
        !           253:   else {
        !           254:     STOQ(e,ne);
        !           255:     if ( !(e%2) ) {
        !           256:       if ( initvp(0,a) ) {
        !           257:         Xmin = 0;
        !           258:         chsgnnum(INF(a),&ai);
        !           259:         if ( compnum(0,ai,SUP(a)) > 0 ) {
        !           260:           Xmax = ai;
        !           261:         } else {
        !           262:           Xmax = SUP(a);
        !           263:         }
        !           264:       } else {
        !           265:         if ( compnum(0,INF(a),0) > 0 ) {
        !           266:           Xmin = INF(a);
        !           267:           Xmax = SUP(a);
        !           268:         } else {
        !           269:           Xmin = SUP(a);
        !           270:           Xmax = INF(a);
        !           271:         }
        !           272:       }
        !           273:     } else {
        !           274:       Xmin = INF(a);
        !           275:       Xmax = SUP(a);
        !           276:     }
        !           277:     if ( ! Xmin )  inf = 0;
        !           278:     else    pwrnum(0,Xmin,ne,&inf);
        !           279:     pwrnum(0,Xmax,ne,&sup);
        !           280:     istoitv(inf,sup,c);
        !           281:   }
1.1       saito     282: }
                    283:
                    284: void chsgnitvp(Itv a, Itv *c)
                    285: {
1.6     ! noro      286:   Num inf,sup;
1.1       saito     287:
1.6     ! noro      288:   if ( !a )
        !           289:     *c = 0;
        !           290:   else if ( NID(a) <= N_B )
        !           291:     chsgnnum(a,c);
        !           292:   else {
        !           293:     chsgnnum(INF((Itv)a),&inf);
        !           294:     chsgnnum(SUP((Itv)a),&sup);
        !           295:     istoitv(inf,sup,c);
        !           296:   }
1.1       saito     297: }
                    298:
                    299: int cmpitvp(Itv a, Itv b)
                    300: {
1.6     ! noro      301:   Num  ai,as,bi,bs;
        !           302:   int  s,t;
1.1       saito     303:
1.6     ! noro      304:   if ( !a ) {
        !           305:     if ( !b || (NID(b)<=N_B) )
        !           306:       return compnum(0,a,b);
        !           307:     else
        !           308:       return -1;
        !           309:   } else if ( !b ) {
        !           310:     if ( !a || (NID(a)<=N_B) )
        !           311:       return compnum(0,a,b);
        !           312:     else
        !           313:       return initvp((Num)b,a);
        !           314:   } else {
        !           315:     itvtois(a,&ai,&as);
        !           316:     itvtois(b,&bi,&bs);
        !           317:     s = compnum(0,ai,bi) ;
        !           318:     t = compnum(0,as,bs) ;
        !           319:     if ( !s && !t ) return 0;
        !           320:     else  return -1;
        !           321:   }
1.1       saito     322: }
                    323:
                    324: void miditvp(Itv a, Num *b)
                    325: {
1.6     ! noro      326:   Num  ai,as;
        !           327:   Num  t;
1.1       saito     328:
1.6     ! noro      329:   if ( ! a ) *b = 0;
        !           330:   else if ( (NID(a) <= N_B) )
        !           331:     *b = (Num)a;
        !           332:   else {
        !           333:     STOQ(2,TWO);
        !           334:     itvtois(a,&ai,&as);
        !           335:     addnum(0,ai,as,&t);
        !           336:     divnum(0,t,TWO,b);
        !           337:   }
1.1       saito     338: }
                    339:
                    340: void cupitvp(Itv a, Itv b, Itv *c)
                    341: {
1.6     ! noro      342:   Num  ai,as,bi,bs;
        !           343:   Num  inf,sup;
1.1       saito     344:
1.6     ! noro      345:   itvtois(a,&ai,&as);
        !           346:   itvtois(b,&bi,&bs);
        !           347:   if ( compnum(0,ai,bi) > 0 )  inf = bi;
        !           348:   else        inf = ai;
        !           349:   if ( compnum(0,as,bs) > 0 )  sup = as;
        !           350:   else        sup = bs;
        !           351:   istoitv(inf,sup,c);
1.1       saito     352: }
                    353:
                    354: void capitvp(Itv a, Itv b, Itv *c)
                    355: {
1.6     ! noro      356:   Num  ai,as,bi,bs;
        !           357:   Num  inf,sup;
1.1       saito     358:
1.6     ! noro      359:   itvtois(a,&ai,&as);
        !           360:   itvtois(b,&bi,&bs);
        !           361:   if ( compnum(0,ai,bi) > 0 )  inf = ai;
        !           362:   else        inf = bi;
        !           363:   if ( compnum(0,as,bs) > 0 )  sup = bs;
        !           364:   else        sup = as;
        !           365:   if ( compnum(0,inf,sup) > 0 )  *c = 0;
        !           366:   else        istoitv(inf,sup,c);
1.1       saito     367: }
                    368:
                    369:
                    370: void widthitvp(Itv a, Num *b)
                    371: {
1.6     ! noro      372:   Num  ai,as;
1.1       saito     373:
1.6     ! noro      374:   if ( ! a ) *b = 0;
        !           375:   else if ( (NID(a) <= N_B) )
        !           376:     *b = (Num)a;
        !           377:   else {
        !           378:     itvtois(a,&ai,&as);
        !           379:     subnum(0,as,ai,b);
        !           380:   }
1.1       saito     381: }
                    382:
                    383: void absitvp(Itv a, Num *b)
                    384: {
1.6     ! noro      385:   Num  ai,as,bi,bs;
1.1       saito     386:
1.6     ! noro      387:   if ( ! a ) *b = 0;
        !           388:   else if ( (NID(a) <= N_B) ) {
        !           389:     if ( compnum(0,a,0) < 0 ) chsgnnum(a,b);
        !           390:     else *b = (Num)a;
        !           391:   } else {
        !           392:     itvtois(a,&ai,&as);
        !           393:     if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
        !           394:     else bi = ai;
        !           395:     if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
        !           396:     else bs = as;
        !           397:     if ( compnum(0,bi,bs) > 0 )  *b = bi;
        !           398:     else        *b = bs;
        !           399:   }
1.1       saito     400: }
                    401:
                    402: void distanceitvp(Itv a, Itv b, Num *c)
                    403: {
1.6     ! noro      404:   Num  ai,as,bi,bs;
        !           405:   Num  di,ds;
        !           406:   Itv  d;
        !           407:
        !           408:   itvtois(a,&ai,&as);
        !           409:   itvtois(b,&bi,&bs);
        !           410:   subnum(0,ai,bi,&di);
        !           411:   subnum(0,as,bs,&ds);
        !           412:   istoitv(di,ds,&d);
        !           413:   absitvp(d,c);
1.1       saito     414: }
                    415:
                    416: #endif

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