version 1.3, 2019/06/04 07:11:23 |
version 1.6, 2019/11/12 10:53:22 |
|
|
/* |
/* |
* $OpenXM: OpenXM_contrib2/asir2018/engine/p-itv.c,v 1.2 2018/09/28 08:20:28 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2018/engine/p-itv.c,v 1.5 2019/11/07 08:50:49 kondoh Exp $ |
*/ |
*/ |
#if defined(INTERVAL) |
#if defined(INTERVAL) |
#include "ca.h" |
#include "ca.h" |
Line 15 Num tobf(Num,int); |
|
Line 15 Num tobf(Num,int); |
|
//int initvp(Num , Itv ); |
//int initvp(Num , Itv ); |
|
|
extern int mpfr_roundmode; |
extern int mpfr_roundmode; |
|
|
extern int zerorewrite; |
extern int zerorewrite; |
|
|
void itvtois(Itv a, Num *inf, Num *sup) |
void itvtois(Itv a, Num *inf, Num *sup) |
Line 24 void itvtois(Itv a, Num *inf, Num *sup) |
|
Line 23 void itvtois(Itv a, Num *inf, Num *sup) |
|
*inf = *sup = (Num)0; |
*inf = *sup = (Num)0; |
else if ( NID(a) <= N_B ) { |
else if ( NID(a) <= N_B ) { |
*inf = (Num)a; *sup = (Num)a; |
*inf = (Num)a; *sup = (Num)a; |
|
} else if ( ITVD(a) ) { |
|
double dinf, dsup; |
|
Real rinf, rsup; |
|
dinf = INF((IntervalDouble)a); dsup = SUP((IntervalDouble)a); |
|
MKReal(dinf, rinf); MKReal(dsup, rsup); |
|
*inf = (Num)rinf; *sup = (Num)rsup; |
} else { |
} else { |
*inf = INF(a); |
*inf = INF(a); |
*sup = SUP(a); |
*sup = SUP(a); |
Line 32 void itvtois(Itv a, Num *inf, Num *sup) |
|
Line 37 void itvtois(Itv a, Num *inf, Num *sup) |
|
|
|
void istoitv(Num inf, Num sup, Itv *rp) |
void istoitv(Num inf, Num sup, Itv *rp) |
{ |
{ |
Num i, s; |
Num ii, ss; |
Num ni, ns; |
Num ni, ns; |
Itv c; |
Itv c; |
int type=0; |
int type=0; |
Line 50 void istoitv(Num inf, Num sup, Itv *rp) |
|
Line 55 void istoitv(Num inf, Num sup, Itv *rp) |
|
|
|
current_roundmode = mpfr_roundmode; |
current_roundmode = mpfr_roundmode; |
if ( !ns ) { |
if ( !ns ) { |
s = 0; |
ss = 0; |
} |
} |
else if ( NID( ns )==N_B ) { |
else if ( NID( ns )==N_B ) { |
type = 1; |
type = 1; |
|
|
mpfr_roundmode = MPFR_RNDU; |
mpfr_roundmode = MPFR_RNDU; |
s = tobf(ns, DEFAULTPREC); |
ss = tobf(ns, DEFAULTPREC); |
//ToBf(sup, (BF *)&s); |
//ToBf(sup, (BF *)&s); |
} else { |
} else { |
s = ns; |
ss = ns; |
} |
} |
if ( !ni ) { |
if ( !ni ) { |
i = 0; |
ii = 0; |
} |
} |
else if ( NID( ni )==N_B ) { |
else if ( NID( ni )==N_B ) { |
type = 1; |
type = 1; |
|
|
mpfr_roundmode = MPFR_RNDD; |
mpfr_roundmode = MPFR_RNDD; |
i = tobf(inf, DEFAULTPREC); |
ii = tobf(ni, DEFAULTPREC); |
//ToBf(inf, (BF *)&i); |
//ToBf(inf, (BF *)&i); |
} else { |
} else { |
i = ni; |
ii = ni; |
} |
} |
mpfr_roundmode = current_roundmode; |
mpfr_roundmode = current_roundmode; |
|
|
Line 82 void istoitv(Num inf, Num sup, Itv *rp) |
|
Line 87 void istoitv(Num inf, Num sup, Itv *rp) |
|
} else { |
} else { |
NEWItvP(c); |
NEWItvP(c); |
} |
} |
INF(c) = i; SUP(c) = s; |
INF(c) = ii; SUP(c) = ss; |
|
|
if (zerorewrite && initvp(0,c) ) { |
*rp = c; |
*rp = 0; |
|
zerorewriteCount++; |
ZEROREWRITE |
} else { |
|
*rp = c; |
|
} |
|
} |
} |
|
|
void additvp(Itv a, Itv b, Itv *c) |
void additvp(Itv a, Itv b, Itv *c) |
Line 272 void divitvp(Itv a, Itv b, Itv *c) |
|
Line 274 void divitvp(Itv a, Itv b, Itv *c) |
|
|
|
void pwritvp(Itv a, Num e, Itv *c) |
void pwritvp(Itv a, Num e, Itv *c) |
{ |
{ |
int ei; |
long ei; |
Itv t; |
Itv t; |
|
|
if ( !e ) |
if ( !e ) |
Line 288 void pwritvp(Itv a, Num e, Itv *c) |
|
Line 290 void pwritvp(Itv a, Num e, Itv *c) |
|
error("pwritv : can't calculate a fractional power"); |
error("pwritv : can't calculate a fractional power"); |
#endif |
#endif |
} else { |
} else { |
ei = QTOS((Q)e); |
//ei = QTOS((Q)e); |
|
ei = mpz_get_si(BDY((Q)e)); |
pwritv0p(a,ei,&t); |
pwritv0p(a,ei,&t); |
if ( SGN((Q)e) < 0 ) |
// if ( SGN((Q)e) < 0 ) |
|
if ( sgnq((Q)e) < 0 ) |
divnum(0,(Num)ONE,(Num)t,(Num *)c); |
divnum(0,(Num)ONE,(Num)t,(Num *)c); |
else |
else |
*c = t; |
*c = t; |
} |
} |
} |
} |
|
|
void pwritv0p(Itv a, int e, Itv *c) |
void pwritv0p(Itv a, long e, Itv *c) |
{ |
{ |
Num inf, sup; |
Num inf, sup; |
Num ai,Xmin,Xmax; |
Num ai,Xmin,Xmax; |
Line 306 void pwritv0p(Itv a, int e, Itv *c) |
|
Line 310 void pwritv0p(Itv a, int e, Itv *c) |
|
if ( e == 1 ) |
if ( e == 1 ) |
*c = a; |
*c = a; |
else { |
else { |
STOQ(e,ne); |
STOZ(e,ne); |
if ( !(e%2) ) { |
if ( !(e%2) ) { |
if ( initvp(0,a) ) { |
if ( initvp(0,a) ) { |
Xmin = 0; |
Xmin = 0; |
Line 394 void miditvp(Itv a, Num *b) |
|
Line 398 void miditvp(Itv a, Num *b) |
|
else if ( (NID(a) <= N_B) ) |
else if ( (NID(a) <= N_B) ) |
*b = (Num)a; |
*b = (Num)a; |
else { |
else { |
STOQ(2,TWO); |
//STOZ(2,TWO); |
itvtois(a,&ai,&as); |
itvtois(a,&ai,&as); |
addnum(0,ai,as,&t); |
addnum(0,ai,as,&t); |
divnum(0,t,(Num)TWO,b); |
divnum(0,t,(Num)TWO,b); |
Line 462 void absitvp(Itv a, Num *b) |
|
Line 466 void absitvp(Itv a, Num *b) |
|
else *b = bs; |
else *b = bs; |
} |
} |
} |
} |
|
|
|
void absintvalp(Itv a, Itv *c) |
|
{ |
|
Num ai,as,mai; |
|
Num inf,sup; |
|
|
|
itvtois(a,&ai,&as); |
|
if ( compnum(0,as,0 ) < 0 ) { |
|
chsgnnum(ai, &sup); |
|
chsgnnum(as, &inf); |
|
} else if ( compnum(0,ai,0) < 0 ) { |
|
inf = 0; |
|
chsgnnum(ai, &mai); |
|
if ( compnum(0,as,mai ) > 0 ) sup = as; |
|
else sup = mai; |
|
} else { |
|
inf = ai; |
|
sup = as; |
|
} |
|
istoitv(inf,sup,c); |
|
} |
|
|
|
|
void distanceitvp(Itv a, Itv b, Num *c) |
void distanceitvp(Itv a, Itv b, Num *c) |
{ |
{ |