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

1.1       saito       1: /*
1.4     ! saito       2:  * $OpenXM: OpenXM_contrib2/asir2000/engine/t-itv.c,v 1.3 2003/02/14 22:29:09 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: {
                     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.1       saito     231:                GEN pa,pe,z;
                    232:                int ltop,lbot;
                    233:
                    234:                ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
                    235:                z = gerepile(ltop,lbot,gpui(pa,pe,prec));
                    236:                patori(z,c); cgiv(z);
                    237: #else
                    238:                error("pwritv : can't calculate a fractional power");
                    239: #endif
                    240:        } else {
                    241:                ei = QTOS((Q)e);
                    242:                pwritv0p(a,ei,&t);
                    243:                if ( SGN((Q)e) < 0 )
                    244:                        divnum(0,(Num)ONE,(Num)t,c);
                    245:                else
                    246:                        *c = t;
                    247:        }
                    248: }
                    249:
                    250: void pwritv0p(Itv a, int e, Itv *c)
                    251: {
                    252:        Num inf, sup;
                    253:        Num ai,Xmin,Xmax;
                    254:        Q       ne;
                    255:
                    256:        if ( e == 1 )
                    257:                *c = a;
                    258:        else {
                    259:                STOQ(e,ne);
                    260:                if ( !(e%2) ) {
                    261:                        if ( initvp(0,a) ) {
                    262:                                Xmin = 0;
                    263:                                chsgnnum(INF(a),&ai);
                    264:                                if ( compnum(0,ai,SUP(a)) > 0 ) {
                    265:                                        Xmax = ai;
                    266:                                } else {
                    267:                                        Xmax = SUP(a);
                    268:                                }
                    269:                        } else {
                    270:                                if ( compnum(0,INF(a),0) > 0 ) {
                    271:                                        Xmin = INF(a);
                    272:                                        Xmax = SUP(a);
                    273:                                } else {
                    274:                                        Xmin = SUP(a);
                    275:                                        Xmax = INF(a);
                    276:                                }
                    277:                        }
                    278:                } else {
                    279:                        Xmin = INF(a);
                    280:                        Xmax = SUP(a);
                    281:                }
                    282:                if ( ! Xmin )   inf = 0;
                    283:                else            pwrnum(0,Xmin,ne,&inf);
                    284:                pwrnum(0,Xmax,ne,&sup);
                    285:                istoitv(inf,sup,c);
                    286:        }
                    287: }
                    288:
                    289: void chsgnitvp(Itv a, Itv *c)
                    290: {
                    291:        Num inf,sup;
                    292:
                    293:        if ( !a )
                    294:                *c = 0;
                    295:        else if ( NID(a) <= N_B )
                    296:                chsgnnum(a,c);
                    297:        else {
                    298:                chsgnnum(INF((Itv)a),&inf);
                    299:                chsgnnum(SUP((Itv)a),&sup);
                    300:                istoitv(inf,sup,c);
                    301:        }
                    302: }
                    303:
                    304: int cmpitvp(Itv a, Itv b)
                    305: {
                    306:        Num     ai,as,bi,bs;
                    307:        int     s,t;
                    308:
                    309:        if ( !a ) {
                    310:                if ( !b || (NID(b)<=N_B) )
                    311:                        return compnum(0,a,b);
                    312:                else
                    313:                        return -1;
                    314:        } else if ( !b ) {
                    315:                if ( !a || (NID(a)<=N_B) )
                    316:                        return compnum(0,a,b);
                    317:                else
                    318:                        return initvp((Num)b,a);
                    319:        } else {
                    320:                itvtois(a,&ai,&as);
                    321:                itvtois(b,&bi,&bs);
                    322:                s = compnum(0,ai,bi) ;
                    323:                t = compnum(0,as,bs) ;
                    324:                if ( !s && !t ) return 0;
                    325:                else  return -1;
                    326:        }
                    327: }
                    328:
                    329: void miditvp(Itv a, Num *b)
                    330: {
                    331:        Num     ai,as;
                    332:        Num     t;
                    333:
                    334:        if ( ! a ) *b = 0;
                    335:        else if ( (NID(a) <= N_B) )
                    336:                *b = (Num)a;
                    337:        else {
                    338:                STOQ(2,TWO);
                    339:                itvtois(a,&ai,&as);
                    340:                addnum(0,ai,as,&t);
                    341:                divnum(0,t,TWO,b);
                    342:        }
                    343: }
                    344:
                    345: void cupitvp(Itv a, Itv b, Itv *c)
                    346: {
                    347:        Num     ai,as,bi,bs;
                    348:        Num     inf,sup;
                    349:
                    350:        itvtois(a,&ai,&as);
                    351:        itvtois(b,&bi,&bs);
                    352:        if ( compnum(0,ai,bi) > 0 )     inf = bi;
                    353:        else                            inf = ai;
                    354:        if ( compnum(0,as,bs) > 0 )     sup = as;
                    355:        else                            sup = bs;
                    356:        istoitv(inf,sup,c);
                    357: }
                    358:
                    359: void capitvp(Itv a, Itv b, Itv *c)
                    360: {
                    361:        Num     ai,as,bi,bs;
                    362:        Num     inf,sup;
                    363:
                    364:        itvtois(a,&ai,&as);
                    365:        itvtois(b,&bi,&bs);
                    366:        if ( compnum(0,ai,bi) > 0 )     inf = ai;
                    367:        else                            inf = bi;
                    368:        if ( compnum(0,as,bs) > 0 )     sup = bs;
                    369:        else                            sup = as;
                    370:        if ( compnum(0,inf,sup) > 0 )   *c = 0;
                    371:        else                            istoitv(inf,sup,c);
                    372: }
                    373:
                    374:
                    375: void widthitvp(Itv a, Num *b)
                    376: {
                    377:        Num     ai,as;
                    378:
                    379:        if ( ! a ) *b = 0;
                    380:        else if ( (NID(a) <= N_B) )
                    381:                *b = (Num)a;
                    382:        else {
                    383:                itvtois(a,&ai,&as);
                    384:                subnum(0,as,ai,b);
                    385:        }
                    386: }
                    387:
                    388: void absitvp(Itv a, Num *b)
                    389: {
                    390:        Num     ai,as,bi,bs;
                    391:
                    392:        if ( ! a ) *b = 0;
                    393:        else if ( (NID(a) <= N_B) ) {
                    394:                if ( compnum(0,a,0) < 0 ) chsgnnum(a,b);
                    395:                else *b = (Num)a;
                    396:        } else {
                    397:                itvtois(a,&ai,&as);
                    398:                if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
                    399:                else bi = ai;
                    400:                if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
                    401:                else bs = as;
                    402:                if ( compnum(0,bi,bs) > 0 )     *b = bi;
                    403:                else                            *b = bs;
                    404:        }
                    405: }
                    406:
                    407: void distanceitvp(Itv a, Itv b, Num *c)
                    408: {
                    409:        Num     ai,as,bi,bs;
                    410:        Num     di,ds;
                    411:        Itv     d;
                    412:
                    413:        itvtois(a,&ai,&as);
                    414:        itvtois(b,&bi,&bs);
                    415:        subnum(0,ai,bi,&di);
                    416:        subnum(0,as,bs,&ds);
                    417:        istoitv(di,ds,&d);
                    418:        absitvp(d,c);
                    419: }
                    420:
                    421: #endif

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