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

Annotation of OpenXM_contrib2/asir2018/engine/p-itv.c, Revision 1.6

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

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