version 1.3, 2019/06/04 07:11:23 |
version 1.5, 2019/11/12 10:53:22 |
|
|
/* |
/* |
* $OpenXM: OpenXM_contrib2/asir2018/engine/f-itv.c,v 1.2 2018/09/28 08:20:28 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2018/engine/f-itv.c,v 1.4 2019/10/17 03:03:12 kondoh Exp $ |
*/ |
*/ |
#if defined(INTERVAL) |
#if defined(INTERVAL) |
#include "ca.h" |
#include "ca.h" |
Line 120 Num tobf(Num,int); |
|
Line 120 Num tobf(Num,int); |
|
|
|
#define BFPREC(a) (((BF)(a))->body->_mpfr_prec) |
#define BFPREC(a) (((BF)(a))->body->_mpfr_prec) |
|
|
|
double mpfr2dblDown(mpfr_t a) |
|
{ |
|
return mpfr_get_d(a,MPFR_RNDD); |
|
} |
|
|
|
double mpfr2dblUp(mpfr_t a) |
|
{ |
|
return mpfr_get_d(a,MPFR_RNDU); |
|
} |
|
|
|
|
|
void toInterval(Num a, int prec, int type, Num *rp) |
|
{ |
|
if ( ! a ) { |
|
*rp = 0; |
|
} else if (type == EvalIntervalDouble) { |
|
if (NID(a)==N_C) { |
|
double inf,sup; |
|
C z; |
|
IntervalDouble re, im; |
|
|
|
if ( ! ((C)a)->r ) { |
|
re = 0; |
|
} else { |
|
inf = toRealDown(((C)a)->r); |
|
sup = toRealUp(((C)a)->r); |
|
MKIntervalDouble(inf,sup,re); |
|
} |
|
if ( ! ((C)a)->i ) { |
|
im = 0; |
|
} else { |
|
inf = toRealDown(((C)a)->i); |
|
sup = toRealUp(((C)a)->i); |
|
MKIntervalDouble(inf,sup,im); |
|
} |
|
if ( !re && !im ) |
|
z = 0; |
|
else |
|
reimtocplx((Num)re,(Num)im,(Num *)&z); |
|
*rp = (Num)z; |
|
} else { |
|
double inf,sup; |
|
IntervalDouble c; |
|
|
|
inf = toRealDown(a); |
|
sup = toRealUp(a); |
|
|
|
MKIntervalDouble(inf,sup,c); |
|
*rp = (Num) c; |
|
} |
|
} else if (type == EvalIntervalBigFloat) { |
|
if (NID(a)==N_C) { |
|
Num ai,as; |
|
Num inf,sup; |
|
C z; |
|
IntervalBigFloat re, im; |
|
int current_roundmode; |
|
|
|
current_roundmode = mpfr_roundmode; |
|
|
|
if ( ! ((C)a)->r ) |
|
re = 0; |
|
else { |
|
itvtois((Itv)((C)a)->r,&ai,&as); |
|
mpfr_roundmode = MPFR_RNDD; |
|
inf = tobf(ai, prec); |
|
mpfr_roundmode = MPFR_RNDU; |
|
sup = tobf(as, prec); |
|
istoitv(inf,sup,(Itv *)&re); |
|
} |
|
|
|
if ( ! ((C)a)->i ) |
|
im = 0; |
|
else { |
|
itvtois((Itv)((C)a)->i,&ai,&as); |
|
mpfr_roundmode = MPFR_RNDD; |
|
inf = tobf(ai, prec); |
|
mpfr_roundmode = MPFR_RNDU; |
|
sup = tobf(as, prec); |
|
istoitv(inf,sup,(Itv *)&im); |
|
} |
|
|
|
mpfr_roundmode = current_roundmode; |
|
reimtocplx((Num)re,(Num)im,(Num *)&z); |
|
*rp = (Num)z; |
|
} else { |
|
Num ai,as; |
|
Num inf,sup; |
|
IntervalBigFloat c; |
|
int current_roundmode; |
|
|
|
itvtois((Itv)a,&ai,&as); |
|
|
|
current_roundmode = mpfr_roundmode; |
|
mpfr_roundmode = MPFR_RNDD; |
|
inf = tobf(ai, prec); |
|
mpfr_roundmode = MPFR_RNDU; |
|
sup = tobf(as, prec); |
|
istoitv(inf,sup,(Itv *)&c); |
|
mpfr_roundmode = current_roundmode; |
|
*rp = (Num) c; |
|
} |
|
} else { |
|
error("toInterval: not supported types."); |
|
*rp = 0; |
|
} |
|
} |
|
|
|
|
void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp) |
void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp) |
{ |
{ |
Num ai,as,bi,bs; |
Num ai,as,bi,bs; |
Line 387 void divitvf(IntervalBigFloat a, IntervalBigFloat b, I |
|
Line 495 void divitvf(IntervalBigFloat a, IntervalBigFloat b, I |
|
|
|
void pwritvf(Itv a, Num e, Itv *c) |
void pwritvf(Itv a, Num e, Itv *c) |
{ |
{ |
int ei; |
long ei; |
Itv t; |
Itv t; |
|
|
if ( !e ) |
if ( !e ) |
Line 403 void pwritvf(Itv a, Num e, Itv *c) |
|
Line 511 void pwritvf(Itv a, Num e, Itv *c) |
|
error("pwritvf() : can't calculate a fractional power"); |
error("pwritvf() : 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)); |
if (ei<0) ei = -ei; |
if (ei<0) ei = -ei; |
pwritv0f(a,ei,&t); |
pwritv0f(a,ei,&t); |
if ( SGN((Q)e) < 0 ) |
// if ( SGN((Q)e) < 0 ) |
|
if ( sgnq((Q)e) < 0 ) |
divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */ |
divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */ |
else |
else |
*c = t; |
*c = t; |
} |
} |
} |
} |
|
|
void pwritv0f(Itv a, int e, Itv *c) |
void pwritv0f(Itv a, long e, Itv *c) |
{ |
{ |
Num inf, sup; |
Num inf, sup; |
Num ai,Xmin,Xmax; |
Num ai,Xmin,Xmax; |
Line 424 void pwritv0f(Itv a, int e, Itv *c) |
|
Line 534 void pwritv0f(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 533 void miditvf(Itv a, Num *b) |
|
Line 643 void miditvf(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); |
//STOQ(2,TWO); |
itvtois(a,&ai,&as); |
itvtois(a,&ai,&as); |
addbf(ai,as,&t); |
addbf(ai,as,&t); |
divbf(t,(Num)TWO,b); |
divbf(t,(Num)TWO,b); |