=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/itvnum.c,v retrieving revision 1.8 retrieving revision 1.13 diff -u -p -r1.8 -r1.13 --- OpenXM_contrib2/asir2000/builtin/itvnum.c 2015/08/08 14:19:41 1.8 +++ OpenXM_contrib2/asir2000/builtin/itvnum.c 2019/11/12 10:52:04 1.13 @@ -1,14 +1,20 @@ /* - * $OpenXM: OpenXM_contrib2/asir2000/builtin/itvnum.c,v 1.7 2014/05/12 16:54:40 saito Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/itvnum.c,v 1.12 2019/06/04 07:11:22 kondoh Exp $ */ #include "ca.h" #include "parse.h" #include "version.h" +#if !defined(ANDROID) #include "../plot/ifplot.h" +#endif -#if defined(INTERVAL) +#include +#include +// in engine/bf.c +Num tobf(Num,int); +#if defined(INTERVAL) static void Pitv(NODE, Obj *); static void Pitvd(NODE, Obj *); static void Pitvbf(NODE, Obj *); @@ -22,12 +28,21 @@ static void Pcup(NODE, Obj *); static void Pcap(NODE, Obj *); static void Pwidth(NODE, Obj *); static void Pdistance(NODE, Obj *); -static void Pitvversion(Q *); -void miditvp(Itv,Num *); -void absitvp(Itv,Num *); -int initvd(Num,IntervalDouble); -int initvp(Num,Itv); -int itvinitvp(Itv,Itv); +static void Pitvversion(NODE, Q *); +static void PzeroRewriteMode(NODE, Obj *); +static void PzeroRewriteCountClear(NODE, Obj *); +static void PzeroRewriteCount(NODE, Obj *); +//void miditvp(Itv,Num *); +//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 static void Pprintmode(NODE, Obj *); @@ -35,471 +50,559 @@ static void Pprintmode(NODE, Obj *); static void ccalc(double **, struct canvas *, int); static void Pifcheck(NODE, Obj *); -#if defined(__osf__) && 0 -int end; +#if defined(__osf__) && 0 +int end; #endif struct ftab interval_tab[] = { - {"printmode",Pprintmode,1}, + {"printmode",Pprintmode,1}, #if defined(INTERVAL) - {"itvd",Pitvd,-2}, - {"intvald",Pitvd,-2}, - {"itv",Pitv,-2}, - {"intval",Pitv,-2}, - {"itvbf",Pitvbf,-2}, - {"intvalbf",Pitvbf,-2}, - {"inf",Pinf,1}, - {"sup",Psup,1}, - {"absintval",Pabsitv,1}, - {"disintval",Pdisjitv,2}, - {"inintval",Pinitv,2}, - {"cup",Pcup,2}, - {"cap",Pcap,2}, - {"mid",Pmid,1}, - {"width",Pwidth,1}, - {"diam",Pwidth,1}, - {"distance",Pdistance,2}, - {"iversion",Pitvversion,0}, + {"itvd",Pitvd,-2}, + {"intvald",Pitvd,-2}, + {"itv",Pitv,-2}, + {"intval",Pitv,-2}, + {"itvbf",Pitvbf,-2}, + {"intvalbf",Pitvbf,-2}, + {"inf",Pinf,1}, + {"sup",Psup,1}, + {"absintval",Pabsitv,1}, + {"disintval",Pdisjitv,2}, + {"inintval",Pinitv,2}, + {"cup",Pcup,2}, + {"cap",Pcap,2}, + {"mid",Pmid,1}, + {"width",Pwidth,1}, + {"diam",Pwidth,1}, + {"distance",Pdistance,2}, + {"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 */ - {"ifcheck",Pifcheck,-7}, + {"ifcheck",Pifcheck,-7}, #endif - {0,0,0}, + {0,0,0}, }; +extern int mpfr_roundmode; + #if defined(INTERVAL) /* plot time check */ static void Pifcheck(NODE arg, Obj *rp) { - Q m2,p2,s_id; - NODE defrange; - LIST xrange,yrange,range[2],list,geom; - VL vl,vl0; - V v[2],av[2]; - int ri,i,j,sign; - P poly; - P var; - NODE n,n0; - Obj t; + Q m2,p2,s_id; + NODE defrange; + LIST xrange,yrange,range[2],list,geom; + VL vl,vl0; + V v[2],av[2]; + int ri,i,j,sign; + P poly; + P var; + NODE n,n0; + Obj t; - struct canvas *can; - MAT m; - pointer **mb; - double **tabe, *px, *px1, *px2; - Q one; - int width, height, ix, iy; - int id; + struct canvas *can; + MAT m; + pointer **mb; + double **tabe, *px, *px1, *px2; + Q one; + int width, height, ix, iy; + int id; - STOQ(-2,m2); STOQ(2,p2); - STOQ(1,one); - MKNODE(n,p2,0); MKNODE(defrange,m2,n); - poly = 0; vl = 0; geom = 0; ri = 0; - v[0] = v[1] = 0; - for ( ; arg; arg = NEXT(arg) ){ - switch ( OID(BDY(arg)) ) { - case O_P: - poly = (P)BDY(arg); - get_vars_recursive((Obj)poly,&vl); - for(vl0=vl,i=0;vl0;vl0=NEXT(vl0)){ - if(vl0->v->attr==(pointer)V_IND){ - if(i>=2){ - error("ifplot : invalid argument"); - } else { - v[i++]=vl0->v; - } - } - } - break; - case O_LIST: - list = (LIST)BDY(arg); - if ( OID(BDY(BDY(list))) == O_P ) - if ( ri > 1 ) - error("ifplot : invalid argument"); - else - range[ri++] = list; - else - geom = list; - break; - default: - error("ifplot : invalid argument"); break; - } - } - if ( !poly ) error("ifplot : invalid argument"); - switch ( ri ) { - case 0: - if ( !v[1] ) error("ifplot : please specify all variables"); - MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n); - MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n); - break; - case 1: - if ( !v[1] ) error("ifplot : please specify all variables"); - av[0] = VR((P)BDY(BDY(range[0]))); - if ( v[0] == av[0] ) { - xrange = range[0]; - MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n); - } else if ( v[1] == av[0] ) { - MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n); - yrange = range[0]; - } else - error("ifplot : invalid argument"); - break; - case 2: - av[0] = VR((P)BDY(BDY(range[0]))); - av[1] = VR((P)BDY(BDY(range[1]))); - if ( ((v[0] == av[0]) && (!v[1] || v[1] == av[1])) || - ((v[0] == av[1]) && (!v[1] || v[1] == av[0])) ) { - xrange = range[0]; yrange = range[1]; - } else error("ifplot : invalid argument"); - break; - default: - error("ifplot : cannot happen"); break; - } - can = canvas[id = search_canvas()]; - if ( !geom ) { - width = 300; - height = 300; - can->width = 300; - can->height = 300; - } else { - can->width = QTOS((Q)BDY(BDY(geom))); - can->height = QTOS((Q)BDY(NEXT(BDY(geom)))); - width = can->width; - height = can->height; - } - if ( xrange ) { - n = BDY(xrange); can->vx = VR((P)BDY(n)); n = NEXT(n); - can->qxmin = (Q)BDY(n); n = NEXT(n); can->qxmax = (Q)BDY(n); - can->xmin = ToReal(can->qxmin); can->xmax = ToReal(can->qxmax); - } - if ( yrange ) { - n = BDY(yrange); can->vy = VR((P)BDY(n)); n = NEXT(n); - can->qymin = (Q)BDY(n); n = NEXT(n); can->qymax = (Q)BDY(n); - can->ymin = ToReal(can->qymin); can->ymax = ToReal(can->qymax); - } - can->wname = "ifcheck"; - can->formula = poly; - tabe = (double **)ALLOCA((width+1)*sizeof(double *)); - for ( i = 0; i <= width; i++ ) - tabe[i] = (double *)ALLOCA((height+1)*sizeof(double)); - for(i=0;i<=width;i++)for(j=0;j<=height;j++)tabe[i][j]=0; - ccalc(tabe,can,0); - MKMAT(m,width,height); - mb = BDY(m); - for( ix=0; ix= 0 ){ - if ( (tabe[ix+1][iy] <= 0) || - (tabe[ix][iy+1] <= 0 ) || - (tabe[ix+1][iy+1] <= 0 ) ) mb[ix][iy] = (Obj)one; - } else { - if ( (tabe[ix+1][iy] >= 0 ) || - ( tabe[ix][iy+1] >= 0 ) || - ( tabe[ix+1][iy+1] >= 0 )) mb[ix][iy] = (Obj)one; - } - } - } - *rp = (Obj)m; + STOQ(-2,m2); STOQ(2,p2); + STOQ(1,one); + MKNODE(n,p2,0); MKNODE(defrange,m2,n); + poly = 0; vl = 0; geom = 0; ri = 0; + v[0] = v[1] = 0; + for ( ; arg; arg = NEXT(arg) ){ + switch ( OID(BDY(arg)) ) { + case O_P: + poly = (P)BDY(arg); + get_vars_recursive((Obj)poly,&vl); + for(vl0=vl,i=0;vl0;vl0=NEXT(vl0)){ + if(vl0->v->attr==(pointer)V_IND){ + if(i>=2){ + error("ifplot : invalid argument"); + } else { + v[i++]=vl0->v; + } + } + } + break; + case O_LIST: + list = (LIST)BDY(arg); + if ( OID(BDY(BDY(list))) == O_P ) + if ( ri > 1 ) + error("ifplot : invalid argument"); + else + range[ri++] = list; + else + geom = list; + break; + default: + error("ifplot : invalid argument"); break; + } + } + if ( !poly ) error("ifplot : invalid argument"); + switch ( ri ) { + case 0: + if ( !v[1] ) error("ifplot : please specify all variables"); + MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n); + MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n); + break; + case 1: + if ( !v[1] ) error("ifplot : please specify all variables"); + av[0] = VR((P)BDY(BDY(range[0]))); + if ( v[0] == av[0] ) { + xrange = range[0]; + MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n); + } else if ( v[1] == av[0] ) { + MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n); + yrange = range[0]; + } else + error("ifplot : invalid argument"); + break; + case 2: + av[0] = VR((P)BDY(BDY(range[0]))); + av[1] = VR((P)BDY(BDY(range[1]))); + if ( ((v[0] == av[0]) && (!v[1] || v[1] == av[1])) || + ((v[0] == av[1]) && (!v[1] || v[1] == av[0])) ) { + xrange = range[0]; yrange = range[1]; + } else error("ifplot : invalid argument"); + break; + default: + error("ifplot : cannot happen"); break; + } + can = canvas[id = search_canvas()]; + if ( !geom ) { + width = 300; + height = 300; + can->width = 300; + can->height = 300; + } else { + can->width = QTOS((Q)BDY(BDY(geom))); + can->height = QTOS((Q)BDY(NEXT(BDY(geom)))); + width = can->width; + height = can->height; + } + if ( xrange ) { + n = BDY(xrange); can->vx = VR((P)BDY(n)); n = NEXT(n); + can->qxmin = (Q)BDY(n); n = NEXT(n); can->qxmax = (Q)BDY(n); + can->xmin = ToReal(can->qxmin); can->xmax = ToReal(can->qxmax); + } + if ( yrange ) { + n = BDY(yrange); can->vy = VR((P)BDY(n)); n = NEXT(n); + can->qymin = (Q)BDY(n); n = NEXT(n); can->qymax = (Q)BDY(n); + can->ymin = ToReal(can->qymin); can->ymax = ToReal(can->qymax); + } + can->wname = "ifcheck"; + can->formula = poly; + tabe = (double **)ALLOCA((width+1)*sizeof(double *)); + for ( i = 0; i <= width; i++ ) + tabe[i] = (double *)ALLOCA((height+1)*sizeof(double)); + for(i=0;i<=width;i++)for(j=0;j<=height;j++)tabe[i][j]=0; + ccalc(tabe,can,0); + MKMAT(m,width,height); + mb = BDY(m); + for( ix=0; ix= 0 ){ + if ( (tabe[ix+1][iy] <= 0) || + (tabe[ix][iy+1] <= 0 ) || + (tabe[ix+1][iy+1] <= 0 ) ) mb[ix][iy] = (Obj)one; + } else { + if ( (tabe[ix+1][iy] >= 0 ) || + ( tabe[ix][iy+1] >= 0 ) || + ( tabe[ix+1][iy+1] >= 0 )) mb[ix][iy] = (Obj)one; + } + } + } + *rp = (Obj)m; } void ccalc(double **tab,struct canvas *can,int nox) { - double x,y,xmin,ymin,xstep,ystep; - int ix,iy; - Real r,rx,ry; - Obj fr,g; - int w,h; - V vx,vy; - Obj t,s; + double x,y,xmin,ymin,xstep,ystep; + int ix,iy; + Real r,rx,ry; + Obj fr,g; + int w,h; + V vx,vy; + Obj t,s; - MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr); - vx = can->vx; - vy = can->vy; - w = can->width; h = can->height; - xmin = can->xmin; xstep = (can->xmax-can->xmin)/w; - ymin = can->ymin; ystep = (can->ymax-can->ymin)/h; - MKReal(1.0,rx); MKReal(1.0,ry); - for( ix = 0, x = xmin; ix < w+1 ; ix++, x += xstep ) { - BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t); - devalr(CO,t,&g); - for( iy = 0, y = ymin; iy < h+1 ; iy++, y += ystep ) { - BDY(ry) = y; - substr(CO,0,g,vy,y?(Obj)ry:0,&t); - devalr(CO,t,&s); - tab[ix][iy] = ToReal(s); - } - } + MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr); + vx = can->vx; + vy = can->vy; + w = can->width; h = can->height; + xmin = can->xmin; xstep = (can->xmax-can->xmin)/w; + ymin = can->ymin; ystep = (can->ymax-can->ymin)/h; + MKReal(1.0,rx); MKReal(1.0,ry); + for( ix = 0, x = xmin; ix < w+1 ; ix++, x += xstep ) { + BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t); + devalr(CO,t,&g); + for( iy = 0, y = ymin; iy < h+1 ; iy++, y += ystep ) { + BDY(ry) = y; + substr(CO,0,g,vy,y?(Obj)ry:0,&t); + devalr(CO,t,&s); + tab[ix][iy] = ToReal(s); + } + } } /* end plot time check */ 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; static void Pitv(NODE arg, Obj *rp) { - Num a, i, s; - Itv c; - double inf, sup; + Num a, i, s; + Itv c; + double inf, sup; #if 1 - if ( bigfloat ) - Pitvbf(arg, rp); - else - Pitvd(arg,rp); + if ( bigfloat ) + Pitvbf(arg, rp); + else + Pitvd(arg,rp); #else - asir_assert(ARG0(arg),O_N,"itv"); - if ( argc(arg) > 1 ) { - asir_assert(ARG1(arg),O_N,"itv"); - istoitv((Num)ARG0(arg),(Num)ARG1(arg),&c); - } else { - a = (Num)ARG0(arg); - if ( ! a ) { - *rp = 0; - return; - } - else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat) { - *rp = (Obj)a; - return; - } - else if ( NID(a) == N_IntervalDouble ) { - inf = INF((IntervalDouble)a); - sup = SUP((IntervalDouble)a); - double2bf(inf, (BF *)&i); - double2bf(sup, (BF *)&s); - istoitv(i,s,&c); - } - else istoitv(a,a,&c); - } - if ( NID( c ) == N_IntervalBigFloat ) - addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp); - else *rp = (Obj)c; + asir_assert(ARG0(arg),O_N,"itv"); + if ( argc(arg) > 1 ) { + asir_assert(ARG1(arg),O_N,"itv"); + istoitv((Num)ARG0(arg),(Num)ARG1(arg),&c); + } else { + a = (Num)ARG0(arg); + if ( ! a ) { + *rp = 0; + return; + } + else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat) { + *rp = (Obj)a; + return; + } + else if ( NID(a) == N_IntervalDouble ) { + inf = INF((IntervalDouble)a); + sup = SUP((IntervalDouble)a); + double2bf(inf, (BF *)&i); + double2bf(sup, (BF *)&s); + istoitv(i,s,&c); + } + else istoitv(a,a,&c); + } + if ( NID( c ) == N_IntervalBigFloat ) + addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp); + else *rp = (Obj)c; #endif } static void Pitvbf(NODE arg, Obj *rp) { - Num a, i, s; - Itv c; - BF ii,ss; - double inf, sup; + Num a, i, s; + IntervalBigFloat c; + Num ii,ss; + Real di, ds; + double inf, sup; + int current_roundmode; - asir_assert(ARG0(arg),O_N,"intvalbf"); - a = (Num)ARG0(arg); - if ( argc(arg) > 1 ) { - asir_assert(ARG1(arg),O_N,"intvalbf"); - i = (Num)ARG0(arg); - s = (Num)ARG1(arg); - ToBf(i, &ii); - ToBf(s, &ss); - istoitv((Num)ii,(Num)ss,&c); - } else { - if ( ! a ) { - *rp = 0; - return; - } - else if ( NID(a) == N_IP ) { - itvtois((Itv)a, &i, &s); - ToBf(i, &ii); - ToBf(s, &ss); - istoitv((Num)ii,(Num)ss,&c); - } - else if ( NID(a) == N_IntervalBigFloat) { - *rp = (Obj)a; - return; - } - else if ( NID(a) == N_IntervalDouble ) { - inf = INF((IntervalDouble)a); - sup = SUP((IntervalDouble)a); - double2bf(inf, (BF *)&i); - double2bf(sup, (BF *)&s); - istoitv(i,s,&c); - } - else { - ToBf(a, (BF *)&i); - istoitv(i,i,&c); - } - } - if ( c && OID( c ) == O_N && NID( c ) == N_IntervalBigFloat ) - addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp); - else *rp = (Obj)c; + asir_assert(ARG0(arg),O_N,"intvalbf"); + a = (Num)ARG0(arg); + if ( argc(arg) > 1 ) { + asir_assert(ARG1(arg),O_N,"intvalbf"); + + i = (Num)ARG0(arg); + s = (Num)ARG1(arg); + current_roundmode = mpfr_roundmode; + mpfr_roundmode = MPFR_RNDD; + 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 { + if ( ! a ) { + *rp = 0; + return; + } + else if ( NID(a) == N_IP ) { + itvtois((Itv)a, &i, &s); + current_roundmode = mpfr_roundmode; + mpfr_roundmode = MPFR_RNDD; + 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) { + *rp = (Obj)a; + return; + } + else if ( NID(a) == N_IntervalDouble ) { + inf = INF((IntervalDouble)a); + sup = SUP((IntervalDouble)a); + current_roundmode = mpfr_roundmode; + //double2bf(inf, (BF *)&i); + //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 { + current_roundmode = mpfr_roundmode; + 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 ) +// addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp); +// else *rp = (Obj)c; + *rp = (Obj)c; } static void Pitvd(NODE arg, Obj *rp) { - double inf, sup; - Num a, a0, a1, t; - Itv ia; - IntervalDouble d; + double inf, sup; + Num a, a0, a1, t; + Itv ia; + IntervalDouble d; - asir_assert(ARG0(arg),O_N,"intvald"); - a0 = (Num)ARG0(arg); - if ( argc(arg) > 1 ) { - asir_assert(ARG1(arg),O_N,"intvald"); - a1 = (Num)ARG1(arg); - } else { - if ( a0 && OID(a0)==O_N && NID(a0)==N_IntervalDouble ) { - inf = INF((IntervalDouble)a0); - sup = SUP((IntervalDouble)a0); - MKIntervalDouble(inf,sup,d); - *rp = (Obj)d; - return; - } - a1 = (Num)ARG0(arg); - } - if ( compnum(0,a0,a1) > 0 ) { - t = a0; a0 = a1; a1 = t; - } - inf = ToRealDown(a0); - sup = ToRealUp(a1); - MKIntervalDouble(inf,sup,d); - *rp = (Obj)d; + asir_assert(ARG0(arg),O_N,"intvald"); + a0 = (Num)ARG0(arg); + if ( argc(arg) > 1 ) { + asir_assert(ARG1(arg),O_N,"intvald"); + a1 = (Num)ARG1(arg); + } else { + if ( a0 && OID(a0)==O_N && NID(a0)==N_IntervalDouble ) { + inf = INF((IntervalDouble)a0); + sup = SUP((IntervalDouble)a0); + MKIntervalDouble(inf,sup,d); + *rp = (Obj)d; + return; + } + a1 = (Num)ARG0(arg); + } + if ( compnum(0,a0,a1) > 0 ) { + t = a0; a0 = a1; a1 = t; + } + inf = toRealDown(a0); + sup = toRealUp(a1); + MKIntervalDouble(inf,sup,d); + *rp = (Obj)d; } static void Pinf(NODE arg, Obj *rp) { - Num a, i, s; - Real r; - double d; + Num a, i, s; + Real r; + double d; - a = (Num)ARG0(arg); - if ( ! a ) { - *rp = 0; - } else if ( OID(a) == O_N ) { - switch ( NID(a) ) { - case N_IntervalDouble: - d = INF((IntervalDouble)a); - MKReal(d, r); - *rp = (Obj)r; - break; - case N_IP: - case N_IntervalBigFloat: - case N_IntervalQuad: - itvtois((Itv)ARG0(arg),&i,&s); - *rp = (Obj)i; - break; - default: - *rp = (Obj)a; - break; - } - } else { - *rp = (Obj)a; - } + a = (Num)ARG0(arg); + if ( ! a ) { + *rp = 0; + } else if ( OID(a) == O_N ) { + switch ( NID(a) ) { + case N_IntervalDouble: + d = INF((IntervalDouble)a); + MKReal(d, r); + *rp = (Obj)r; + break; + case N_IP: + case N_IntervalBigFloat: + case N_IntervalQuad: + itvtois((Itv)ARG0(arg),&i,&s); + *rp = (Obj)i; + break; + default: + *rp = (Obj)a; + break; + } + } else { + *rp = (Obj)a; + } } static void Psup(NODE arg, Obj *rp) { - Num a, i, s; - Real r; - double d; + Num a, i, s; + Real r; + double d; - a = (Num)ARG0(arg); - if ( ! a ) { - *rp = 0; - } else if ( OID(a) == O_N ) { - switch ( NID(a) ) { - case N_IntervalDouble: - d = SUP((IntervalDouble)a); - MKReal(d, r); - *rp = (Obj)r; - break; - case N_IP: - case N_IntervalBigFloat: - case N_IntervalQuad: - itvtois((Itv)ARG0(arg),&i,&s); - *rp = (Obj)s; - break; - default: - *rp = (Obj)a; - break; - } - } else { - *rp = (Obj)a; - } + a = (Num)ARG0(arg); + if ( ! a ) { + *rp = 0; + } else if ( OID(a) == O_N ) { + switch ( NID(a) ) { + case N_IntervalDouble: + d = SUP((IntervalDouble)a); + MKReal(d, r); + *rp = (Obj)r; + break; + case N_IP: + case N_IntervalBigFloat: + case N_IntervalQuad: + itvtois((Itv)ARG0(arg),&i,&s); + *rp = (Obj)s; + break; + default: + *rp = (Obj)a; + break; + } + } else { + *rp = (Obj)a; + } } static void Pmid(NODE arg, Obj *rp) { - Num a, s; - Real r; - double d; + Num a, s; + Real r; + double d; - a = (Num)ARG0(arg); - if ( ! a ) { - *rp = 0; - } else switch (OID(a)) { - case O_N: - if ( NID(a) == N_IntervalDouble ) { - d = ( INF((IntervalDouble)a)+SUP((IntervalDouble)a) ) / 2.0; - MKReal(d, r); - *rp = (Obj)r; - } else if ( NID(a) == N_IntervalQuad ) { - error("mid: not supported operation"); - *rp = 0; - } else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat ) { - miditvp((Itv)ARG0(arg),&s); - *rp = (Obj)s; - } else { - *rp = (Obj)a; - } - break; + a = (Num)ARG0(arg); + if ( ! a ) { + *rp = 0; + } else switch (OID(a)) { + case O_N: + if ( NID(a) == N_IntervalDouble ) { + d = ( INF((IntervalDouble)a)+SUP((IntervalDouble)a) ) / 2.0; + MKReal(d, r); + *rp = (Obj)r; + } else if ( NID(a) == N_IntervalQuad ) { + error("mid: not supported operation"); + *rp = 0; + } else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat ) { + miditvp((Itv)ARG0(arg),&s); + *rp = (Obj)s; + } else { + *rp = (Obj)a; + } + break; #if 0 - case O_P: - case O_R: - case O_LIST: - case O_VECT: - case O_MAT: + case O_P: + case O_R: + case O_LIST: + case O_VECT: + case O_MAT: #endif - default: - *rp = (Obj)a; - break; - } + default: + *rp = (Obj)a; + break; + } } static void Pcup(NODE arg, Obj *rp) { - Itv s; - Num a, b; + Itv s; + Num a, b; - asir_assert(ARG0(arg),O_N,"cup"); - asir_assert(ARG1(arg),O_N,"cup"); - a = (Num)ARG0(arg); - b = (Num)ARG1(arg); - if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) { - cupitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp); - } else { - cupitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s); - *rp = (Obj)s; - } + asir_assert(ARG0(arg),O_N,"cup"); + asir_assert(ARG1(arg),O_N,"cup"); + a = (Num)ARG0(arg); + b = (Num)ARG1(arg); + if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) { + cupitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp); + } else { + cupitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s); + *rp = (Obj)s; + } } static void Pcap(NODE arg, Obj *rp) { - Itv s; - Num a, b; + Itv s; + Num a, b; - asir_assert(ARG0(arg),O_N,"cap"); - asir_assert(ARG1(arg),O_N,"cap"); - a = (Num)ARG0(arg); - b = (Num)ARG1(arg); - if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) { - capitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp); - } else { - capitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s); - *rp = (Obj)s; - } + asir_assert(ARG0(arg),O_N,"cap"); + asir_assert(ARG1(arg),O_N,"cap"); + a = (Num)ARG0(arg); + b = (Num)ARG1(arg); + if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) { + capitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp); + } else { + capitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s); + *rp = (Obj)s; + } } static void @@ -507,19 +610,19 @@ Pwidth(arg,rp) NODE arg; Obj *rp; { - Num s; - Num a; + Num s; + Num a; - asir_assert(ARG0(arg),O_N,"width"); - a = (Num)ARG0(arg); - if ( ! a ) { - *rp = 0; - } else if ( NID(a) == N_IntervalDouble ) { - widthitvd((IntervalDouble)a, (Num *)rp); - } else { - widthitvp((Itv)ARG0(arg),&s); - *rp = (Obj)s; - } + asir_assert(ARG0(arg),O_N,"width"); + a = (Num)ARG0(arg); + if ( ! a ) { + *rp = 0; + } else if ( NID(a) == N_IntervalDouble ) { + widthitvd((IntervalDouble)a, (Num *)rp); + } else { + widthitvp((Itv)ARG0(arg),&s); + *rp = (Obj)s; + } } static void @@ -527,19 +630,19 @@ Pabsitv(arg,rp) NODE arg; Obj *rp; { - Num s; - Num a, b; + Num s; + Num a, b; - asir_assert(ARG0(arg),O_N,"absitv"); - a = (Num)ARG0(arg); - if ( ! a ) { - *rp = 0; - } else if ( NID(a) == N_IntervalDouble ) { - absitvd((IntervalDouble)a, (Num *)rp); - } else { - absitvp((Itv)ARG0(arg),&s); - *rp = (Obj)s; - } + asir_assert(ARG0(arg),O_N,"absitv"); + a = (Num)ARG0(arg); + if ( ! a ) { + *rp = 0; + } else if ( NID(a) == N_IntervalDouble ) { + absitvd((IntervalDouble)a, (Num *)rp); + } else { + absitvp((Itv)ARG0(arg),&s); + *rp = (Obj)s; + } } static void @@ -547,19 +650,19 @@ Pdistance(arg,rp) NODE arg; Obj *rp; { - Num s; - Num a, b; + Num s; + Num a, b; - asir_assert(ARG0(arg),O_N,"distance"); - asir_assert(ARG1(arg),O_N,"distance"); - a = (Num)ARG0(arg); - b = (Num)ARG1(arg); - if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) { - distanceitvd((IntervalDouble)a, (IntervalDouble)b, (Num *)rp); - } else { - distanceitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s); - *rp = (Obj)s; - } + asir_assert(ARG0(arg),O_N,"distance"); + asir_assert(ARG1(arg),O_N,"distance"); + a = (Num)ARG0(arg); + b = (Num)ARG1(arg); + if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) { + distanceitvd((IntervalDouble)a, (IntervalDouble)b, (Num *)rp); + } else { + distanceitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s); + *rp = (Obj)s; + } } static void @@ -567,30 +670,30 @@ Pinitv(arg,rp) NODE arg; Obj *rp; { - int s; - Q q; + int s; + Q q; - asir_assert(ARG0(arg),O_N,"intval"); - asir_assert(ARG1(arg),O_N,"intval"); - if ( ! ARG1(arg) ) { - if ( ! ARG0(arg) ) s = 1; - else s = 0; - } - else if ( NID(ARG1(arg)) == N_IntervalDouble ) { - s = initvd((Num)ARG0(arg),(IntervalDouble)ARG1(arg)); + asir_assert(ARG0(arg),O_N,"intval"); + asir_assert(ARG1(arg),O_N,"intval"); + if ( ! ARG1(arg) ) { + if ( ! ARG0(arg) ) s = 1; + else s = 0; + } + else if ( NID(ARG1(arg)) == N_IntervalDouble ) { + s = initvd((Num)ARG0(arg),(IntervalDouble)ARG1(arg)); - } else if ( NID(ARG1(arg)) == N_IP || NID(ARG1(arg)) == N_IntervalBigFloat ) { - if ( ! ARG0(arg) ) s = initvp((Num)ARG0(arg),(Itv)ARG1(arg)); - else if ( NID(ARG0(arg)) == N_IP ) { - s = itvinitvp((Itv)ARG0(arg),(Itv)ARG1(arg)); - } else { - s = initvp((Num)ARG0(arg),(Itv)ARG1(arg)); - } - } else { - s = ! compnum(0,(Num)ARG0(arg),(Num)ARG1(arg)); - } - STOQ(s,q); - *rp = (Obj)q; + } else if ( NID(ARG1(arg)) == N_IP || NID(ARG1(arg)) == N_IntervalBigFloat ) { + if ( ! ARG0(arg) ) s = initvp((Num)ARG0(arg),(Itv)ARG1(arg)); + else if ( NID(ARG0(arg)) == N_IP ) { + s = itvinitvp((Itv)ARG0(arg),(Itv)ARG1(arg)); + } else { + s = initvp((Num)ARG0(arg),(Itv)ARG1(arg)); + } + } else { + s = ! compnum(0,(Num)ARG0(arg),(Num)ARG1(arg)); + } + STOQ(s,q); + *rp = (Obj)q; } static void @@ -598,64 +701,791 @@ Pdisjitv(arg,rp) NODE arg; Obj *rp; { - Itv s; + Itv s; - asir_assert(ARG0(arg),O_N,"disjitv"); - asir_assert(ARG1(arg),O_N,"disjitv"); - error("disjitv: not implemented yet"); - if ( ! s ) *rp = 0; - else *rp = (Obj)ONE; + asir_assert(ARG0(arg),O_N,"disjitv"); + asir_assert(ARG1(arg),O_N,"disjitv"); + error("disjitv: not implemented yet"); + if ( ! s ) *rp = 0; + 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 -extern int printmode; +extern int printmode; -static void pprintmode( void ) +static void pprintmode( void ) { - switch (printmode) { + switch (printmode) { #if defined(INTERVAL) - case MID_PRINTF_E: - fprintf(stderr,"Interval printing mode is a mitpoint type.\n"); + case MID_PRINTF_E: + fprintf(stderr,"Interval printing mode is a mitpoint type.\n"); #endif - case PRINTF_E: - fprintf(stderr,"Printf's double printing mode is \"%%.16e\".\n"); - break; + case PRINTF_E: + fprintf(stderr,"Printf's double printing mode is \"%%.16e\".\n"); + break; #if defined(INTERVAL) - case MID_PRINTF_G: - fprintf(stderr,"Interval printing mode is a mitpoint type.\n"); + case MID_PRINTF_G: + fprintf(stderr,"Interval printing mode is a mitpoint type.\n"); #endif - default: - case PRINTF_G: - fprintf(stderr,"Printf's double printing mode is \"%%g\".\n"); - break; - } -#if defined(__MINGW32__) || defined(__MINGW64__) - fflush(stderr); -#endif + default: + case PRINTF_G: + fprintf(stderr,"Printf's double printing mode is \"%%g\".\n"); + break; + } } static void Pprintmode(NODE arg, Obj *rp) { - int l; - Q a, r; + int l; + Q a, r; - a = (Q)ARG0(arg); - if(!a||(NUM(a)&&INT(a))){ - l=QTOS(a); - if ( l < 0 ) l = 0; + a = (Q)ARG0(arg); + if(!a||(NUM(a)&&INT(a))){ + l=QTOS(a); + if ( l < 0 ) l = 0; #if defined(INTERVAL) - else if ( l > MID_PRINTF_E ) l = 0; + else if ( l > MID_PRINTF_E ) l = 0; #else - else if ( l > PRINTF_E ) l = 0; + else if ( l > PRINTF_E ) l = 0; #endif - STOQ(printmode,r); - *rp = (Obj)r; - printmode = l; - pprintmode(); + STOQ(printmode,r); + *rp = (Obj)r; + printmode = l; + pprintmode(); + } else { + *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 { - *rp = 0; + (*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