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

Annotation of OpenXM_contrib2/asir2000/engine/p-itv.c, Revision 1.9

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

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