[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.5

1.1       saito       1: /*
1.5     ! ohara       2:  * $OpenXM: OpenXM_contrib2/asir2000/engine/t-itv.c,v 1.4 2003/10/20 07:18:42 saito 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: {
                     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) )
1.2       kondoh     53:                additvd((Num)a,(Num)b,(IntervalDouble *)c);
1.1       saito      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) )
1.2       kondoh     76:                subitvd((Num)a,(Num)b,(IntervalDouble *)c);
1.1       saito      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) )
1.2       kondoh     97:                mulitvd((Num)a,(Num)b,(IntervalDouble *)c);
1.1       saito      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) )
1.2       kondoh    179:                divitvd((Num)a,(Num)b,(IntervalDouble *)c);
1.1       saito     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) ) {
1.3       ohara     230: #if defined(PARI) && 0
1.5     ! ohara     231:                gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1       saito     232: #else
                    233:                error("pwritv : can't calculate a fractional power");
                    234: #endif
                    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:        }
                    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 {
                    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:        }
                    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 {
                    333:                STOQ(2,TWO);
                    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>