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

Annotation of OpenXM_contrib2/asir2018/engine/t-itv.c, Revision 1.2

1.1       noro        1: /*
1.2     ! noro        2:  * $OpenXM: OpenXM_contrib2/asir2018/engine/t-itv.c,v 1.1 2018/09/19 05:45:07 noro Exp $
1.1       noro        3: */
                      4: #if defined(INTERVAL)
                      5: #include "ca.h"
                      6: #include "base.h"
                      7: #if defined(PARI)
                      8: #include "genpari.h"
                      9: #endif
                     10:
                     11: void itvtois(Itv a, Num *inf, Num *sup)
                     12: {
                     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:   }
                     21: }
                     22:
                     23: void istoitv(Num inf, Num sup, Itv *rp)
                     24: {
                     25:   Itv c;
                     26:
                     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:   }
                     38: }
                     39:
                     40: void additvp(Itv a, Itv b, Itv *c)
                     41: {
                     42:   Num  ai,as,bi,bs;
                     43:   Num  inf,sup;
                     44:
                     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:   }
                     61: }
                     62:
                     63: void subitvp(Itv a, Itv b, Itv *c)
                     64: {
                     65:   Num  ai,as,bi,bs;
                     66:   Num  inf,sup;
                     67:
                     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:   }
                     84: }
                     85:
                     86: void mulitvp(Itv a, Itv b, Itv *c)
                     87: {
                     88:   Num  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t;
                     89:   int  ba,bb;
                     90:
                     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:   }
                    144: }
                    145:
                    146: int     initvp(Num a, Itv b)
                    147: {
                    148:   Num inf, sup;
                    149:
                    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;
                    154: }
                    155:
                    156: int     itvinitvp(Itv a, Itv b)
                    157: {
                    158:   Num ai,as,bi,bs;
                    159:
                    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;
                    164: }
                    165:
                    166: void divitvp(Itv a, Itv b, Itv *c)
                    167: {
                    168:   Num  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t;
                    169:   int  ba,bb;
                    170:
                    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:   }
                    216: }
                    217:
                    218: void pwritvp(Itv a, Num e, Itv *c)
                    219: {
                    220:   int ei;
                    221:   Itv t;
                    222:
                    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) ) {
                    230: #if defined(PARI) && 0
                    231:     gpui_ri((Obj)a,(Obj)c,(Obj *)c);
                    232: #else
                    233:     error("pwritv : can't calculate a fractional power");
                    234: #endif
                    235:   } else {
1.2     ! noro      236:     ei = ZTOS((Q)e);
1.1       noro      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:   }
                    243: }
                    244:
                    245: void pwritv0p(Itv a, int e, Itv *c)
                    246: {
                    247:   Num inf, sup;
                    248:   Num ai,Xmin,Xmax;
                    249:   Q  ne;
                    250:
                    251:   if ( e == 1 )
                    252:     *c = a;
                    253:   else {
1.2     ! noro      254:     STOZ(e,ne);
1.1       noro      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:   }
                    282: }
                    283:
                    284: void chsgnitvp(Itv a, Itv *c)
                    285: {
                    286:   Num inf,sup;
                    287:
                    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:   }
                    297: }
                    298:
                    299: int cmpitvp(Itv a, Itv b)
                    300: {
                    301:   Num  ai,as,bi,bs;
                    302:   int  s,t;
                    303:
                    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:   }
                    322: }
                    323:
                    324: void miditvp(Itv a, Num *b)
                    325: {
                    326:   Num  ai,as;
                    327:   Num  t;
                    328:
                    329:   if ( ! a ) *b = 0;
                    330:   else if ( (NID(a) <= N_B) )
                    331:     *b = (Num)a;
                    332:   else {
1.2     ! noro      333:     STOZ(2,TWO);
1.1       noro      334:     itvtois(a,&ai,&as);
                    335:     addnum(0,ai,as,&t);
                    336:     divnum(0,t,TWO,b);
                    337:   }
                    338: }
                    339:
                    340: void cupitvp(Itv a, Itv b, Itv *c)
                    341: {
                    342:   Num  ai,as,bi,bs;
                    343:   Num  inf,sup;
                    344:
                    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);
                    352: }
                    353:
                    354: void capitvp(Itv a, Itv b, Itv *c)
                    355: {
                    356:   Num  ai,as,bi,bs;
                    357:   Num  inf,sup;
                    358:
                    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);
                    367: }
                    368:
                    369:
                    370: void widthitvp(Itv a, Num *b)
                    371: {
                    372:   Num  ai,as;
                    373:
                    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:   }
                    381: }
                    382:
                    383: void absitvp(Itv a, Num *b)
                    384: {
                    385:   Num  ai,as,bi,bs;
                    386:
                    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:   }
                    400: }
                    401:
                    402: void distanceitvp(Itv a, Itv b, Num *c)
                    403: {
                    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);
                    414: }
                    415:
                    416: #endif

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