version 1.11, 2018/03/29 01:32:50 |
version 1.13, 2019/11/12 10:52:04 |
|
|
/* |
/* |
* $OpenXM: OpenXM_contrib2/asir2000/builtin/itvnum.c,v 1.10 2016/06/29 08:16:11 ohara Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/builtin/itvnum.c,v 1.12 2019/06/04 07:11:22 kondoh Exp $ |
*/ |
*/ |
|
|
#include "ca.h" |
#include "ca.h" |
|
|
#include "../plot/ifplot.h" |
#include "../plot/ifplot.h" |
#endif |
#endif |
|
|
#if defined(INTERVAL) |
#include <mpfr.h> |
|
#include <mpfi.h> |
|
// in engine/bf.c |
|
Num tobf(Num,int); |
|
|
|
#if defined(INTERVAL) |
static void Pitv(NODE, Obj *); |
static void Pitv(NODE, Obj *); |
static void Pitvd(NODE, Obj *); |
static void Pitvd(NODE, Obj *); |
static void Pitvbf(NODE, Obj *); |
static void Pitvbf(NODE, Obj *); |
Line 24 static void Pcup(NODE, Obj *); |
|
Line 28 static void Pcup(NODE, Obj *); |
|
static void Pcap(NODE, Obj *); |
static void Pcap(NODE, Obj *); |
static void Pwidth(NODE, Obj *); |
static void Pwidth(NODE, Obj *); |
static void Pdistance(NODE, Obj *); |
static void Pdistance(NODE, Obj *); |
static void Pitvversion(Q *); |
static void Pitvversion(NODE, Q *); |
void miditvp(Itv,Num *); |
static void PzeroRewriteMode(NODE, Obj *); |
void absitvp(Itv,Num *); |
static void PzeroRewriteCountClear(NODE, Obj *); |
int initvd(Num,IntervalDouble); |
static void PzeroRewriteCount(NODE, Obj *); |
int initvp(Num,Itv); |
//void miditvp(Itv,Num *); |
int itvinitvp(Itv,Itv); |
//void absitvp(Itv,Num *); |
|
//int initvd(Num,IntervalDouble); |
|
//int initvp(Num,Itv); |
|
//int itvinitvp(Itv,Itv); |
|
static void Pevalitv(NODE, Obj *); |
|
static void Pevalitvd(NODE, Obj *); |
|
void Ppi_itvd(NODE, Obj *); |
|
void Pe_itvd(NODE, Obj *); |
|
void Psinitv(NODE, Obj *); |
|
void Psinitvd(NODE, Obj *); |
#endif |
#endif |
static void Pprintmode(NODE, Obj *); |
static void Pprintmode(NODE, Obj *); |
|
|
Line 61 struct ftab interval_tab[] = { |
|
Line 74 struct ftab interval_tab[] = { |
|
{"width",Pwidth,1}, |
{"width",Pwidth,1}, |
{"diam",Pwidth,1}, |
{"diam",Pwidth,1}, |
{"distance",Pdistance,2}, |
{"distance",Pdistance,2}, |
{"iversion",Pitvversion,0}, |
{"iversion",Pitvversion,-1}, |
|
{"intvalversion",Pitvversion,-1}, |
|
{"zerorewritemode",PzeroRewriteMode,-1}, |
|
{"zeroRewriteMode",PzeroRewriteMode,-1}, |
|
{"zeroRewriteCountClear",PzeroRewriteCountClear,-1}, |
|
{"zeroRewriteCount",PzeroRewriteCount,-1}, |
|
/* eval */ |
|
{"evalitv", Pevalitv, -2}, |
|
{"evalitvd", Pevalitvd, 1}, |
|
/* math */ |
|
{"piitvd", Pitvbf_pi, -1}, |
|
{"eitvd", Pitvbf_e, -1}, |
|
|
|
{"piitv", Pitvbf_pi, -1}, |
|
{"eitv", Pitvbf_e, -1}, |
|
#if 0 |
|
{"factorialitv",Pfactorialitv,1}, |
|
{"factorialitvd",Pfactorialitvd,1}, |
|
#endif |
|
|
|
{"absitv", Pitvbf_abs, -2}, |
|
{"absitvd", Pitvbf_abs, -2}, |
|
|
|
{"logitv", Pitvbf_log, -2}, |
|
{"logitvd", Pitvbf_log, -2}, |
|
{"expitv", Pitvbf_exp, -2}, |
|
{"expitvd", Pitvbf_exp, -2}, |
|
{"powitv", Pitvbf_pow, -3}, |
|
{"powitvd", Pitvbf_pow, -3}, |
|
|
|
{"sinitv", Pitvbf_sin, -2}, |
|
{"sinitvd", Pitvd_sin, -2}, |
|
|
|
{"cositv", Pitvbf_cos, -2}, |
|
{"cositvd", Pitvd_cos, -2}, |
|
{"tanitv", Pitvbf_tan, -2}, |
|
{"tanitvd", Pitvd_tan, -2}, |
|
{"asinitv", Pitvbf_asin, -2}, |
|
{"asinitvd", Pitvd_asin, -2}, |
|
{"acositv", Pitvbf_acos, -2}, |
|
{"acositvd", Pitvd_acos, -2}, |
|
{"atanitv", Pitvbf_atan, -2}, |
|
{"atanitvd", Pitvd_atan, -2}, |
|
{"sinhitv", Pitvbf_sinh, -2}, |
|
{"sinhitvd", Pitvd_sinh, -2}, |
|
{"coshitv", Pitvbf_cosh, -2}, |
|
{"coshitvd", Pitvd_cosh, -2}, |
|
{"tanhitv", Pitvbf_tanh, -2}, |
|
{"tanhitvd", Pitvd_tanh, -2}, |
|
{"asinhitv", Pitvbf_asinh, -2}, |
|
{"asinhitvd", Pitvd_asinh, -2}, |
|
{"acoshitv", Pitvbf_acosh, -2}, |
|
{"acoshitvd", Pitvd_acosh, -2}, |
|
{"atanhitv", Pitvbf_atanh, -2}, |
|
{"atanhitvd", Pitvd_atanh, -2}, |
|
|
/* plot time check */ |
/* plot time check */ |
{"ifcheck",Pifcheck,-7}, |
{"ifcheck",Pifcheck,-7}, |
#endif |
#endif |
{0,0,0}, |
{0,0,0}, |
}; |
}; |
|
|
|
extern int mpfr_roundmode; |
|
|
#if defined(INTERVAL) |
#if defined(INTERVAL) |
|
|
/* plot time check */ |
/* plot time check */ |
Line 235 void ccalc(double **tab,struct canvas *can,int nox) |
|
Line 305 void ccalc(double **tab,struct canvas *can,int nox) |
|
/* end plot time check */ |
/* end plot time check */ |
|
|
static void |
static void |
Pitvversion(Q *rp) |
Pitvversion(NODE arg, Q *rp) |
{ |
{ |
STOQ(ASIR_VERSION, *rp); |
STOQ(INT_ASIR_VERSION, *rp); |
} |
} |
|
|
extern int bigfloat; |
extern int bigfloat; |
|
|
Pitvbf(NODE arg, Obj *rp) |
Pitvbf(NODE arg, Obj *rp) |
{ |
{ |
Num a, i, s; |
Num a, i, s; |
Itv c; |
IntervalBigFloat c; |
BF ii,ss; |
Num ii,ss; |
|
Real di, ds; |
double inf, sup; |
double inf, sup; |
|
int current_roundmode; |
|
|
asir_assert(ARG0(arg),O_N,"intvalbf"); |
asir_assert(ARG0(arg),O_N,"intvalbf"); |
a = (Num)ARG0(arg); |
a = (Num)ARG0(arg); |
if ( argc(arg) > 1 ) { |
if ( argc(arg) > 1 ) { |
asir_assert(ARG1(arg),O_N,"intvalbf"); |
asir_assert(ARG1(arg),O_N,"intvalbf"); |
|
|
i = (Num)ARG0(arg); |
i = (Num)ARG0(arg); |
s = (Num)ARG1(arg); |
s = (Num)ARG1(arg); |
ToBf(i, &ii); |
current_roundmode = mpfr_roundmode; |
ToBf(s, &ss); |
mpfr_roundmode = MPFR_RNDD; |
istoitv((Num)ii,(Num)ss,&c); |
ii = tobf(i, DEFAULTPREC); |
|
mpfr_roundmode = MPFR_RNDU; |
|
ss = tobf(s, DEFAULTPREC); |
|
istoitv(ii,ss,(Itv *)&c); |
|
// MKIntervalBigFloat((BF)ii,(BF)ss,c); |
|
// ToBf(s, &ss); |
|
mpfr_roundmode = current_roundmode; |
} else { |
} else { |
if ( ! a ) { |
if ( ! a ) { |
*rp = 0; |
*rp = 0; |
Line 308 Pitvbf(NODE arg, Obj *rp) |
|
Line 387 Pitvbf(NODE arg, Obj *rp) |
|
} |
} |
else if ( NID(a) == N_IP ) { |
else if ( NID(a) == N_IP ) { |
itvtois((Itv)a, &i, &s); |
itvtois((Itv)a, &i, &s); |
ToBf(i, &ii); |
current_roundmode = mpfr_roundmode; |
ToBf(s, &ss); |
mpfr_roundmode = MPFR_RNDD; |
istoitv((Num)ii,(Num)ss,&c); |
ii = tobf(i, DEFAULTPREC); |
|
mpfr_roundmode = MPFR_RNDU; |
|
ss = tobf(s, DEFAULTPREC); |
|
istoitv(ii,ss,(Itv *)&c); |
|
// MKIntervalBigFloat((BF)ii,(BF)ss,c); |
|
mpfr_roundmode = current_roundmode; |
} |
} |
else if ( NID(a) == N_IntervalBigFloat) { |
else if ( NID(a) == N_IntervalBigFloat) { |
*rp = (Obj)a; |
*rp = (Obj)a; |
Line 319 Pitvbf(NODE arg, Obj *rp) |
|
Line 403 Pitvbf(NODE arg, Obj *rp) |
|
else if ( NID(a) == N_IntervalDouble ) { |
else if ( NID(a) == N_IntervalDouble ) { |
inf = INF((IntervalDouble)a); |
inf = INF((IntervalDouble)a); |
sup = SUP((IntervalDouble)a); |
sup = SUP((IntervalDouble)a); |
double2bf(inf, (BF *)&i); |
current_roundmode = mpfr_roundmode; |
double2bf(sup, (BF *)&s); |
//double2bf(inf, (BF *)&i); |
istoitv(i,s,&c); |
//double2bf(sup, (BF *)&s); |
|
mpfr_roundmode = MPFR_RNDD; |
|
MKReal(inf,di); |
|
ii = tobf((Num)di, DEFAULTPREC); |
|
mpfr_roundmode = MPFR_RNDU; |
|
MKReal(sup,ds); |
|
ss = tobf((Num)ds, DEFAULTPREC); |
|
istoitv(ii,ss,(Itv *)&c); |
|
// MKIntervalBigFloat((BF)ii,(BF)ss,c); |
|
mpfr_roundmode = current_roundmode; |
} |
} |
else { |
else { |
ToBf(a, (BF *)&i); |
current_roundmode = mpfr_roundmode; |
istoitv(i,i,&c); |
mpfr_roundmode = MPFR_RNDD; |
|
ii = tobf(a, DEFAULTPREC); |
|
mpfr_roundmode = MPFR_RNDU; |
|
ss = tobf(a, DEFAULTPREC); |
|
//ToBf(a, (BF *)&i); |
|
istoitv(ii,ss,(Itv *)&c); |
|
// MKIntervalBigFloat((BF)ii,(BF)ss,c); |
|
mpfr_roundmode = current_roundmode; |
} |
} |
} |
} |
if ( c && OID( c ) == O_N && NID( c ) == N_IntervalBigFloat ) |
// if ( c && OID( c ) == O_N && NID( c ) == N_IntervalBigFloat ) |
addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp); |
// addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp); |
else *rp = (Obj)c; |
// else *rp = (Obj)c; |
|
*rp = (Obj)c; |
} |
} |
|
|
static void |
static void |
Line 359 Pitvd(NODE arg, Obj *rp) |
|
Line 460 Pitvd(NODE arg, Obj *rp) |
|
if ( compnum(0,a0,a1) > 0 ) { |
if ( compnum(0,a0,a1) > 0 ) { |
t = a0; a0 = a1; a1 = t; |
t = a0; a0 = a1; a1 = t; |
} |
} |
inf = ToRealDown(a0); |
inf = toRealDown(a0); |
sup = ToRealUp(a1); |
sup = toRealUp(a1); |
MKIntervalDouble(inf,sup,d); |
MKIntervalDouble(inf,sup,d); |
*rp = (Obj)d; |
*rp = (Obj)d; |
} |
} |
|
|
else *rp = (Obj)ONE; |
else *rp = (Obj)ONE; |
} |
} |
|
|
|
static void |
|
PzeroRewriteMode(NODE arg, Obj *rp) |
|
{ |
|
Q a, r; |
|
|
|
STOQ(zerorewrite,r); |
|
*rp = (Obj)r; |
|
|
|
if (arg) { |
|
a = (Q)ARG0(arg); |
|
if(!a) { |
|
zerorewrite = 0; |
|
} else if ( (NUM(a)&&INT(a)) ){ |
|
zerorewrite = 1; |
|
} |
|
} |
|
} |
|
|
|
static void |
|
PzeroRewriteCountClear(NODE arg, Obj *rp) |
|
{ |
|
Q a, r; |
|
|
|
STOQ(zerorewriteCount,r); |
|
*rp = (Obj)r; |
|
|
|
if (arg) { |
|
a = (Q)ARG0(arg); |
|
if(a &&(NUM(a)&&INT(a))){ |
|
zerorewriteCount = 0; |
|
} |
|
} |
|
} |
|
|
|
static void |
|
PzeroRewriteCount(NODE arg, Obj *rp) |
|
{ |
|
Q r; |
|
|
|
STOQ(zerorewriteCount,r); |
|
*rp = (Obj)r; |
|
} |
|
|
|
|
#endif |
#endif |
extern int printmode; |
extern int printmode; |
|
|
Line 656 Pprintmode(NODE arg, Obj *rp) |
|
Line 801 Pprintmode(NODE arg, Obj *rp) |
|
*rp = 0; |
*rp = 0; |
} |
} |
} |
} |
|
|
|
#if defined(INTERVAL) |
|
void |
|
Ppi_itvd(NODE arg, Obj *rp) |
|
{ |
|
double inf, sup; |
|
IntervalDouble c; |
|
FPMINUSINF |
|
sscanf("3.1415926535897932384626433832795028841971693993751", "%lf", &inf); |
|
FPPLUSINF |
|
sscanf("3.1415926535897932384626433832795028841971693993752", "%lf", &sup); |
|
FPNEAREST |
|
MKIntervalDouble(inf,sup,c); |
|
*rp = (Obj)c; |
|
} |
|
void |
|
Pe_itvd(NODE arg, Obj *rp) |
|
{ |
|
double inf, sup; |
|
IntervalDouble c; |
|
FPMINUSINF |
|
sscanf( "2.7182818284590452353602874713526624977572470936999", "%lf", &inf); |
|
FPPLUSINF |
|
sscanf( "2.7182818284590452353602874713526624977572470937000", "%lf", &sup); |
|
FPNEAREST |
|
MKIntervalDouble(inf,sup,c); |
|
*rp = (Obj)c; |
|
} |
|
void |
|
Pln2_itvd(NODE arg, Obj *rp) |
|
{ |
|
double inf, sup; |
|
IntervalDouble c; |
|
FPMINUSINF |
|
sscanf( "0.69314718055994530941723212145817656807550013436025", "%lf", &inf); |
|
FPPLUSINF |
|
sscanf( "0.69314718055994530941723212145817656807550013436026", "%lf", &sup); |
|
FPNEAREST |
|
MKIntervalDouble(inf,sup,c); |
|
*rp = (Obj)c; |
|
} |
|
|
|
void mpfi_func(NODE arg, int (*mpfi_f)(), int prec, Obj *rp) |
|
{ |
|
Num a, ii, ss; |
|
Itv c; |
|
BF inf, sup; |
|
int arg1prec; |
|
mpfi_t mpitv, rv; |
|
|
|
|
|
/* |
|
if ( argc(arg) == 2 ) { |
|
prec = QTOS((Q)ARG1(arg))*3.32193; |
|
if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN; |
|
else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX; |
|
} else { |
|
prec = 0; |
|
prec = mpfr_get_default_prec(); |
|
} |
|
*/ |
|
if ( prec > 0 ) arg1prec = prec; |
|
else arg1prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec(); |
|
a = ARG0(arg); |
|
itvtois((Itv)a, &ii, &ss); |
|
|
|
inf = (BF)tobf(ii, arg1prec); |
|
sup = (BF)tobf(ss, arg1prec); |
|
|
|
mpfi_init2(rv,arg1prec); |
|
mpfi_init2(mpitv,arg1prec); |
|
mpfr_set(&(mpitv->left), BDY(inf), MPFR_RNDD); |
|
mpfr_set(&(mpitv->right), BDY(sup), MPFR_RNDU); |
|
|
|
(*mpfi_f)(rv, mpitv); |
|
//mpfi_sin(rv, mpitv); |
|
|
|
MPFRTOBF(&(rv->left), inf); |
|
MPFRTOBF(&(rv->right), sup); |
|
|
|
if ( !cmpbf((Num)inf,0) ) inf = 0; |
|
if ( !cmpbf((Num)sup,0) ) sup = 0; |
|
|
|
if ( inf || sup ) { |
|
istoitv((Num)inf, (Num)sup, &c); |
|
} else { |
|
c = 0; |
|
} |
|
*rp = (Obj)c; |
|
//mpfi_clear(rv); |
|
mpfi_clear(mpitv); |
|
} |
|
|
|
void mpfi_func_d(NODE arg, int (*mpfi_f)(), Obj *rp) |
|
{ |
|
Obj bfv; |
|
Num ii, ss; |
|
IntervalDouble c; |
|
double inf, sup; |
|
|
|
mpfi_func(arg, mpfi_f, 53, &bfv); |
|
itvtois((Itv)bfv, &ii, &ss); |
|
inf = toRealDown(ii); |
|
sup = toRealUp(ss); |
|
MKIntervalDouble(inf,sup,c); |
|
*rp = (Obj)c; |
|
} |
|
|
|
|
|
void |
|
Psinitvd(NODE arg, Obj *rp) |
|
{ |
|
Obj bfv; |
|
Num ii,ss; |
|
IntervalDouble c; |
|
double ai,as,mas, bi,bs; |
|
double inf,sup; |
|
|
|
#if 1 |
|
mpfi_func(arg, mpfi_sin, 53, &bfv); |
|
itvtois((Itv)bfv, &ii, &ss); |
|
inf = toRealDown(ii); |
|
sup = toRealUp(ss); |
|
MKIntervalDouble(inf,sup,c); |
|
*rp = (Obj)c; |
|
#else |
|
a = ARG0(arg); |
|
Num2double(a,&ai,&as); |
|
FPMINUSINF |
|
inf = sin(ai); |
|
FPPLUSINF |
|
sup = sin(as); |
|
FPNEAREST |
|
MKIntervalDouble(inf,sup,c); |
|
*rp = (Obj)c; |
|
#endif |
|
} |
|
|
|
void |
|
Psinitv(NODE arg, Obj *rp) |
|
{ |
|
//Num a; |
|
//Itv c; |
|
//BF inf, sup; |
|
//int prec; |
|
//BF r,re,im; |
|
//mpfi_t mpitv, rv; |
|
|
|
#if 1 |
|
mpfi_func(arg, mpfi_sin, 0, rp); |
|
#else |
|
prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec(); |
|
a = ARG0(arg); |
|
itvtois((Itv)a, (Num *)&inf, (Num *)&sup); |
|
|
|
mpfi_init2(rv,prec); |
|
mpfi_init2(mpitv,prec); |
|
mpfr_set(&(mpitv->left), inf->body, MPFR_RNDD); |
|
mpfr_set(&(mpitv->right), BDY(sup), MPFR_RNDU); |
|
|
|
//(*mpfi_f)(rv, mpitv); |
|
mpfi_sin(rv, mpitv); |
|
|
|
MPFRTOBF(&(rv->left), inf); |
|
MPFRTOBF(&(rv->right), sup); |
|
istoitv((Num)inf, (Num)sup, &c); |
|
*rp = (Obj)c; |
|
mpfi_clear(rv); |
|
mpfi_clear(mpitv); |
|
#endif |
|
} |
|
|
|
|
|
|
|
|
|
//void evalitvr(VL, Obj, int, int, Obj *); |
|
|
|
static void |
|
Pevalitv(NODE arg, Obj *rp) |
|
{ |
|
int prec; |
|
|
|
asir_assert(ARG0(arg),O_R,"evalitv"); |
|
if ( argc(arg) == 2 ) { |
|
long int mpfr_prec_max = MPFR_PREC_MAX; |
|
prec = QTOS((Q)ARG1(arg))*3.32193; |
|
if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN; |
|
//else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX; |
|
else if ( prec > mpfr_prec_max ) prec = mpfr_prec_max; |
|
} else |
|
prec = 0; |
|
evalitvr(CO,(Obj)ARG0(arg),prec,EvalIntervalBigFloat,rp); |
|
} |
|
|
|
static void |
|
Pevalitvd(NODE arg, Obj *rp) |
|
{ |
|
asir_assert(ARG0(arg),O_R,"evalitvd"); |
|
evalitvr(CO,(Obj)ARG0(arg),53,EvalIntervalDouble,rp); |
|
} |
|
|
|
// in parse/puref.c |
|
void instoobj(PFINS ins,Obj *rp); |
|
|
|
// in this |
|
void evalitvr(VL, Obj, int, int, Obj *); |
|
void evalitvp(VL, P, int, int, P *); |
|
void evalitvv(VL, V, int, int, Obj *); |
|
void evalitvins(PFINS, int, int, Obj *); |
|
|
|
|
|
|
|
void evalitvr(VL vl,Obj a,int prec, int type, Obj *c) |
|
{ |
|
Obj nm,dn; |
|
|
|
if ( !a ) |
|
*c = 0; |
|
else { |
|
switch ( OID(a) ) { |
|
case O_N: |
|
toInterval((Num)a, prec, type, (Num *)c); |
|
break; |
|
case O_P: |
|
evalitvp(vl,(P)a,prec,type,(P *)c); |
|
break; |
|
case O_R: |
|
evalitvp(vl,NM((R)a),prec,type,(P *)&nm); |
|
evalitvp(vl,DN((R)a),prec,type,(P *)&dn); |
|
divr(vl,nm,dn,c); |
|
break; |
|
default: |
|
error("evalr : not implemented"); break; |
|
} |
|
} |
|
} |
|
|
|
void evalitvp(VL vl,P p,int prec, int type, P *pr) |
|
{ |
|
P t; |
|
DCP dc,dcr0,dcr; |
|
Obj u; |
|
|
|
if ( !p || NUM(p) ) { |
|
toInterval((Num)p, prec, type, (Num *)pr); |
|
//*pr = p; |
|
} else { |
|
for ( dcr0 = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) { |
|
evalitvp(vl,COEF(dc),prec,type, &t); |
|
if ( t ) { |
|
NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t; |
|
} |
|
} |
|
if ( !dcr0 ) { |
|
*pr = 0; return; |
|
} else { |
|
NEXT(dcr) = 0; MKP(VR(p),dcr0,t); |
|
} |
|
if ( NUM(t) ) { |
|
//*pr = t; |
|
toInterval((Num)t, prec, type, (Num *)pr); |
|
return; |
|
} else if ( (VR(t) != VR(p)) || ((vid)VR(p)->attr != V_PF) ) { |
|
*pr = t; return; |
|
} else { |
|
evalitvv(vl,VR(p),prec,type,&u); |
|
substr(vl,1,(Obj)t,VR(p),u, (Obj *)&t); |
|
if ( t && NUM(t) ) { |
|
toInterval((Num)t, prec, type, (Num *)pr); |
|
} else { |
|
*pr = t; |
|
} |
|
} |
|
} |
|
} |
|
|
|
void evalitvv(VL vl,V v,int prec, int type, Obj *rp) |
|
{ |
|
PFINS ins,tins; |
|
PFAD ad,tad; |
|
PF pf; |
|
P t; |
|
int i; |
|
|
|
if ( (vid)v->attr != V_PF ) { |
|
MKV(v,t); *rp = (Obj)t; |
|
} else { |
|
ins = (PFINS)v->priv; ad = ins->ad; pf = ins->pf; |
|
tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD)); |
|
tins->pf = pf; |
|
for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) { |
|
tad[i].d = ad[i].d; |
|
evalitvr(vl,ad[i].arg,prec,type,&tad[i].arg); |
|
} |
|
evalitvins(tins,prec,type,rp); |
|
} |
|
} |
|
|
|
void evalitvins(PFINS ins,int prec, int type, Obj *rp) |
|
{ |
|
PF pf; |
|
PFINS tins; |
|
PFAD ad,tad; |
|
int i; |
|
Q q; |
|
V v; |
|
P x; |
|
NODE n0,n; |
|
|
|
|
|
pf = ins->pf; ad = ins->ad; |
|
tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD)); |
|
tins->pf = pf; tad = tins->ad; |
|
for ( i = 0; i < pf->argc; i++ ) { |
|
//tad[i].d = ad[i].d; evalr(CO,ad[i].arg,prec,&tad[i].arg); |
|
tad[i].d = ad[i].d; evalitvr(CO,ad[i].arg,prec,type,&tad[i].arg); |
|
} |
|
for ( i = 0; i < pf->argc; i++ ) |
|
if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break; |
|
if ( (i != pf->argc) || !pf->intervalfunc[type] ) { /////////////////////////// |
|
instoobj(tins,rp); |
|
} else { |
|
int IsCPLX = 0; |
|
for ( n0 = 0, i = 0; i < pf->argc; i++ ) { |
|
NEXTNODE(n0,n); BDY(n) = (pointer)tad[i].arg; |
|
if (tad[i].arg && NID(tad[i].arg) == N_C) IsCPLX = 1; |
|
} |
|
if ( prec ) { |
|
NEXTNODE(n0,n); STOQ(prec,q); BDY(n) = (pointer)q; |
|
} |
|
if ( n0 ) NEXT(n) = 0; |
|
|
|
|
|
if ( IsCPLX ) { |
|
instoobj(tins,rp); |
|
} else { |
|
(*pf->intervalfunc[type])(n0,rp); |
|
} |
|
} |
|
} |
|
|
|
|
|
void Pitvbf_pi(NODE arg, Obj *rp) |
|
{ |
|
BF inf, sup; |
|
IntervalBigFloat c; |
|
mpfi_t rv; |
|
int prec; |
|
|
|
prec = (arg) ? QTOS((Q)ARG0(arg)) : 0; //mpfr_get_default_prec(); |
|
prec ? mpfi_init2(rv,prec) : mpfi_init(rv); |
|
|
|
mpfi_const_pi(rv); |
|
|
|
MPFRTOBF(&(rv->left), inf); |
|
MPFRTOBF(&(rv->right), sup); |
|
|
|
MKIntervalBigFloat(inf,sup,c); |
|
|
|
*rp = (Obj)c; |
|
} |
|
|
|
void Pitvd_pi(NODE arg, Obj *rp) |
|
{ |
|
BF bfinf, bfsup; |
|
IntervalDouble c; |
|
mpfi_t rv; |
|
double inf, sup; |
|
|
|
mpfi_init2(rv,53); |
|
|
|
mpfi_const_pi(rv); |
|
|
|
MPFRTOBF(&(rv->left), bfinf); |
|
MPFRTOBF(&(rv->right), bfsup); |
|
|
|
inf = toRealDown((Num)bfinf); |
|
sup = toRealUp((Num)bfsup); |
|
|
|
MKIntervalDouble(inf,sup,c); |
|
|
|
*rp = (Obj)c; |
|
} |
|
|
|
void Pitvbf_e(NODE arg,Obj *rp) |
|
{ |
|
BF inf, sup; |
|
IntervalBigFloat c; |
|
mpfi_t rv; |
|
mpfi_t one; |
|
int prec; |
|
|
|
prec = (arg) ? QTOS((Q)ARG0(arg)) : 0; //mpfr_get_default_prec(); |
|
prec ? mpfi_init2(rv,prec) : mpfi_init(rv); |
|
|
|
mpfi_init(one); |
|
mpfi_set_ui(one,1); |
|
mpfi_exp(rv,one); |
|
|
|
MPFRTOBF(&(rv->left), inf); |
|
MPFRTOBF(&(rv->right), sup); |
|
|
|
MKIntervalBigFloat(inf,sup,c); |
|
//istoitv((Num)inf, (Num)sup, &c); |
|
*rp = (Obj)c; |
|
mpfi_clear(one); |
|
} |
|
|
|
void Pitvd_e(NODE arg, Obj *rp) |
|
{ |
|
BF bfinf, bfsup; |
|
IntervalDouble c; |
|
mpfi_t rv; |
|
mpfi_t one; |
|
double inf, sup; |
|
|
|
mpfi_init2(rv,53); |
|
|
|
mpfi_init2(one, 53); |
|
mpfi_set_ui(one,1); |
|
mpfi_exp(rv,one); |
|
|
|
|
|
MPFRTOBF(&(rv->left), bfinf); |
|
MPFRTOBF(&(rv->right), bfsup); |
|
|
|
inf = toRealDown((Num)bfinf); |
|
sup = toRealUp((Num)bfsup); |
|
|
|
MKIntervalDouble(inf,sup,c); |
|
|
|
*rp = (Obj)c; |
|
} |
|
|
|
void (*pi_itv_ft[])() = {Pitvd_pi, 0, Pitvbf_pi}; |
|
void (*e_itv_ft[])() = {Pitvd_e, 0, Pitvbf_e}; |
|
//void (*sin_itv_ft[])() = {Psinitvd, 0, Psinitv}; |
|
void (*sin_itv_ft[])() = {Pitvd_sin, 0, Pitvbf_sin}; |
|
void (*cos_itv_ft[])() = {Pitvbf_cos, 0, Pitvbf_cos}; |
|
void (*tan_itv_ft[])() = {Pitvd_tan, 0, Pitvbf_tan}; |
|
void (*asin_itv_ft[])() = {Pitvd_asin, 0, Pitvbf_asin}; |
|
void (*acos_itv_ft[])() = {Pitvd_acos, 0, Pitvbf_acos}; |
|
void (*atan_itv_ft[])() = {Pitvd_atan, 0, Pitvbf_atan}; |
|
void (*sinh_itv_ft[])() = {Pitvd_sinh, 0, Pitvbf_sinh}; |
|
void (*cosh_itv_ft[])() = {Pitvd_cosh, 0, Pitvbf_cosh}; |
|
void (*tanh_itv_ft[])() = {Pitvd_tanh, 0, Pitvbf_tanh}; |
|
void (*asinh_itv_ft[])() = {Pitvd_asinh, 0, Pitvbf_asinh}; |
|
void (*acosh_itv_ft[])() = {Pitvd_acosh, 0, Pitvbf_acosh}; |
|
void (*atanh_itv_ft[])() = {Pitvd_atanh, 0, Pitvbf_atanh}; |
|
void (*exp_itv_ft[])() = {Pitvd_exp, 0, Pitvbf_exp}; |
|
void (*log_itv_ft[])() = {Pitvd_log, 0, Pitvbf_log}; |
|
void (*abs_itv_ft[])() = {0}; |
|
void (*pow_itv_ft[])() = {Pitvbf_pow, 0, Pitvbf_pow}; |
|
//void (*pow_itv_ft[])() = {0, 0, 0}; |
|
|
|
|
|
void Pitvbf_sin(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_sin, 0, rp); |
|
} |
|
|
|
void Pitvbf_cos(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_cos, 0, rp); |
|
} |
|
|
|
void Pitvbf_tan(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_tan, 0, rp); |
|
} |
|
|
|
void Pitvbf_asin(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_asin, 0, rp); |
|
} |
|
|
|
void Pitvbf_acos(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_acos, 0, rp); |
|
} |
|
|
|
void Pitvbf_atan(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_atan, 0, rp); |
|
} |
|
|
|
void Pitvbf_sinh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_sinh, 0, rp); |
|
} |
|
|
|
void Pitvbf_cosh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_cosh, 0, rp); |
|
} |
|
|
|
void Pitvbf_tanh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_tanh, 0, rp); |
|
} |
|
|
|
void Pitvbf_asinh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_asinh, 0, rp); |
|
} |
|
|
|
void Pitvbf_acosh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_acosh, 0, rp); |
|
} |
|
|
|
void Pitvbf_atanh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_atanh, 0, rp); |
|
} |
|
|
|
void Pitvbf_exp(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_exp, 0, rp); |
|
} |
|
|
|
void Pitvbf_log(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_log, 0, rp); |
|
} |
|
|
|
void Pitvbf_abs(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func(arg, mpfi_abs, 0, rp); |
|
} |
|
|
|
void Pitvd_sin(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_sin, rp); |
|
} |
|
|
|
void Pitvd_cos(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_cos, rp); |
|
} |
|
|
|
void Pitvd_tan(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_tan, rp); |
|
} |
|
|
|
void Pitvd_asin(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_asin, rp); |
|
} |
|
|
|
void Pitvd_acos(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_acos, rp); |
|
} |
|
|
|
void Pitvd_atan(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_atan, rp); |
|
} |
|
|
|
void Pitvd_sinh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_sinh, rp); |
|
} |
|
|
|
void Pitvd_cosh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_cosh, rp); |
|
} |
|
|
|
void Pitvd_tanh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_tanh, rp); |
|
} |
|
|
|
void Pitvd_asinh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_asinh, rp); |
|
} |
|
|
|
void Pitvd_acosh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_acosh, rp); |
|
} |
|
|
|
void Pitvd_atanh(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_atanh, rp); |
|
} |
|
|
|
void Pitvd_exp(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_exp, rp); |
|
} |
|
|
|
void Pitvd_log(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_log, rp); |
|
} |
|
|
|
void Pitvd_abs(NODE arg,Obj *rp) |
|
{ |
|
mpfi_func_d(arg, mpfi_abs, rp); |
|
} |
|
|
|
/* |
|
void mp_factorial(NODE arg,Num *rp) |
|
{ |
|
struct oNODE arg0; |
|
Num a,a1; |
|
|
|
a = (Num)ARG0(arg); |
|
if ( !a ) *rp = (Num)ONE; |
|
else if ( INT(a) ) Pfac(arg,rp); |
|
else { |
|
addnum(0,a,(Num)ONE,&a1); |
|
arg0.body = (pointer)a1; |
|
arg0.next = arg->next; |
|
Pmpfr_gamma(&arg0,rp); |
|
} |
|
} |
|
*/ |
|
|
|
void Pitvbf_pow(NODE arg,Num *rp) |
|
{ |
|
Num a,e; |
|
BF r,re,im; |
|
C c; |
|
mpc_t mpc,a1,e1; |
|
int prec; |
|
|
|
prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : mpfr_get_default_prec(); |
|
a = ARG0(arg); |
|
e = ARG1(arg); |
|
if ( !e ) { |
|
*rp = (Num)ONE; |
|
} else if ( !a ) { |
|
*rp = 0; |
|
} else if ( NID(a) == N_C || NID(e) == N_C ) { |
|
error("itv_pow() : not support args"); |
|
*rp = 0; |
|
} else { |
|
Num ii, ss; |
|
Itv c; |
|
BF inf, sup; |
|
mpfi_t a_val, e_val, rv; |
|
|
|
mpfi_init2(rv,prec); |
|
|
|
itvtois((Itv)a, &ii, &ss); |
|
inf = (BF)tobf(ii, prec); |
|
sup = (BF)tobf(ss, prec); |
|
mpfi_init2(a_val,prec); |
|
mpfr_set(&(a_val->left), BDY(inf), MPFR_RNDD); |
|
mpfr_set(&(a_val->right), BDY(sup), MPFR_RNDU); |
|
|
|
itvtois((Itv)e, &ii, &ss); |
|
inf = (BF)tobf(ii, prec); |
|
sup = (BF)tobf(ss, prec); |
|
mpfi_init2(e_val,prec); |
|
mpfr_set(&(e_val->left), BDY(inf), MPFR_RNDD); |
|
mpfr_set(&(e_val->right), BDY(sup), MPFR_RNDU); |
|
|
|
|
|
mpfi_log(rv, a_val); |
|
mpfi_mul(a_val, rv, e_val); |
|
mpfi_exp(rv, a_val); |
|
|
|
MPFRTOBF(&(rv->left), inf); |
|
MPFRTOBF(&(rv->right), sup); |
|
|
|
if ( !cmpbf((Num)inf,0) ) inf = 0; |
|
if ( !cmpbf((Num)sup,0) ) sup = 0; |
|
|
|
if ( inf || sup ) { |
|
istoitv((Num)inf, (Num)sup, &c); |
|
} else { |
|
c = 0; |
|
} |
|
*rp = (Num)c; |
|
//mpfi_clear(rv); |
|
mpfi_clear(a_val); |
|
mpfi_clear(e_val); |
|
} |
|
} |
|
|
|
#endif |