Annotation of OpenXM_contrib2/asir2018/builtin/itvnum.c, Revision 1.6
1.1 noro 1: /*
1.6 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2018/builtin/itvnum.c,v 1.5 2019/11/12 10:53:22 kondoh Exp $
1.1 noro 3: */
4:
5: #include "ca.h"
6: #include "parse.h"
7: #include "version.h"
8: #if !defined(ANDROID)
9: #include "../plot/ifplot.h"
10: #endif
11:
1.5 kondoh 12: #include <mpfr.h>
13: #include <mpfi.h>
1.3 kondoh 14: // in engine/bf.c
15: Num tobf(Num,int);
16:
1.1 noro 17: #if defined(INTERVAL)
18: static void Pitv(NODE, Obj *);
19: static void Pitvd(NODE, Obj *);
20: static void Pitvbf(NODE, Obj *);
21: static void Pinf(NODE, Obj *);
22: static void Psup(NODE, Obj *);
23: static void Pmid(NODE, Obj *);
24: static void Pabsitv(NODE, Obj *);
25: static void Pdisjitv(NODE, Obj *);
26: static void Pinitv(NODE, Obj *);
27: static void Pcup(NODE, Obj *);
28: static void Pcap(NODE, Obj *);
29: static void Pwidth(NODE, Obj *);
30: static void Pdistance(NODE, Obj *);
1.3 kondoh 31: static void Pitvversion(NODE, Q *);
32: static void PzeroRewriteMode(NODE, Obj *);
33: static void PzeroRewriteCountClear(NODE, Obj *);
34: static void PzeroRewriteCount(NODE, Obj *);
35: //void miditvp(Itv,Num *);
36: //void absitvp(Itv,Num *);
37: //int initvd(Num,IntervalDouble);
38: //int initvp(Num,Itv);
39: //int itvinitvp(Itv,Itv);
1.5 kondoh 40: static void Pevalitv(NODE, Obj *);
1.6 ! kondoh 41: static void Pevalitvbf(NODE, Obj *);
1.5 kondoh 42: static void Pevalitvd(NODE, Obj *);
1.6 ! kondoh 43:
! 44: static void Pitvbf_pi(NODE ,Obj *);
! 45: static void Pitvbf_e(NODE ,Obj *);
! 46: static void Pitvbf_sin(NODE ,Obj *);
! 47: static void Pitvbf_cos(NODE ,Obj *);
! 48: static void Pitvbf_tan(NODE ,Obj *);
! 49: static void Pitvbf_asin(NODE ,Obj *);
! 50: static void Pitvbf_acos(NODE ,Obj *);
! 51: static void Pitvbf_atan(NODE ,Obj *);
! 52: static void Pitvbf_sinh(NODE ,Obj *);
! 53: static void Pitvbf_cosh(NODE ,Obj *);
! 54: static void Pitvbf_tanh(NODE ,Obj *);
! 55: static void Pitvbf_asinh(NODE ,Obj *);
! 56: static void Pitvbf_acosh(NODE ,Obj *);
! 57: static void Pitvbf_atanh(NODE ,Obj *);
! 58: static void Pitvbf_exp(NODE ,Obj *);
! 59: static void Pitvbf_log(NODE ,Obj *);
! 60: static void Pitvbf_abs(NODE ,Obj *);
! 61: static void Pitvbf_pow(NODE ,Num *);
! 62:
! 63: static void Pitvd_pi(NODE ,Obj *);
! 64: static void Pitvd_e(NODE ,Obj *);
! 65: static void Pitvd_sin(NODE ,Obj *);
! 66: static void Pitvd_cos(NODE ,Obj *);
! 67: static void Pitvd_tan(NODE ,Obj *);
! 68: static void Pitvd_asin(NODE ,Obj *);
! 69: static void Pitvd_acos(NODE ,Obj *);
! 70: static void Pitvd_atan(NODE ,Obj *);
! 71: static void Pitvd_sinh(NODE ,Obj *);
! 72: static void Pitvd_cosh(NODE ,Obj *);
! 73: static void Pitvd_tanh(NODE ,Obj *);
! 74: static void Pitvd_asinh(NODE ,Obj *);
! 75: static void Pitvd_acosh(NODE ,Obj *);
! 76: static void Pitvd_atanh(NODE ,Obj *);
! 77: static void Pitvd_exp(NODE ,Obj *);
! 78: static void Pitvd_log(NODE ,Obj *);
! 79: static void Pitvd_abs(NODE ,Obj *);
! 80: static void Pitvd_pow(NODE ,Num *);
! 81:
! 82: static void Pitv_pi(NODE ,Obj *);
! 83: static void Pitv_e(NODE ,Obj *);
! 84: static void Pitv_sin(NODE ,Obj *);
! 85: static void Pitv_cos(NODE ,Obj *);
! 86: static void Pitv_tan(NODE ,Obj *);
! 87: static void Pitv_asin(NODE ,Obj *);
! 88: static void Pitv_acos(NODE ,Obj *);
! 89: static void Pitv_atan(NODE ,Obj *);
! 90: static void Pitv_sinh(NODE ,Obj *);
! 91: static void Pitv_cosh(NODE ,Obj *);
! 92: static void Pitv_tanh(NODE ,Obj *);
! 93: static void Pitv_asinh(NODE ,Obj *);
! 94: static void Pitv_acosh(NODE ,Obj *);
! 95: static void Pitv_atanh(NODE ,Obj *);
! 96: static void Pitv_exp(NODE ,Obj *);
! 97: static void Pitv_log(NODE ,Obj *);
! 98: static void Pitv_abs(NODE ,Obj *);
! 99: static void Pitv_pow(NODE ,Num *);
1.1 noro 100: #endif
101: static void Pprintmode(NODE, Obj *);
102:
103: /* plot time check func */
104: static void ccalc(double **, struct canvas *, int);
105: static void Pifcheck(NODE, Obj *);
106:
107: #if defined(__osf__) && 0
108: int end;
109: #endif
110:
111: struct ftab interval_tab[] = {
112: {"printmode",Pprintmode,1},
113: #if defined(INTERVAL)
114: {"itvd",Pitvd,-2},
115: {"intvald",Pitvd,-2},
116: {"itv",Pitv,-2},
117: {"intval",Pitv,-2},
118: {"itvbf",Pitvbf,-2},
119: {"intvalbf",Pitvbf,-2},
120: {"inf",Pinf,1},
121: {"sup",Psup,1},
122: {"absintval",Pabsitv,1},
1.6 ! kondoh 123: {"absitv",Pabsitv,1},
1.1 noro 124: {"disintval",Pdisjitv,2},
125: {"inintval",Pinitv,2},
126: {"cup",Pcup,2},
127: {"cap",Pcap,2},
128: {"mid",Pmid,1},
129: {"width",Pwidth,1},
130: {"diam",Pwidth,1},
131: {"distance",Pdistance,2},
1.3 kondoh 132: {"iversion",Pitvversion,-1},
133: {"intvalversion",Pitvversion,-1},
134: {"zerorewritemode",PzeroRewriteMode,-1},
135: {"zeroRewriteMode",PzeroRewriteMode,-1},
136: {"zeroRewriteCountClear",PzeroRewriteCountClear,-1},
137: {"zeroRewriteCount",PzeroRewriteCount,-1},
1.5 kondoh 138: /* eval */
139: {"evalitv", Pevalitv, -2},
1.6 ! kondoh 140: {"evalitvbf", Pevalitvbf, -2},
1.5 kondoh 141: {"evalitvd", Pevalitvd, 1},
142: /* math */
143:
1.6 ! kondoh 144: {"piitv", Pitv_pi, -1},
! 145: {"piitvbf", Pitvbf_pi, -1},
! 146: {"piitvd", Pitvd_pi, -1},
! 147: {"eitv", Pitv_e, -1},
! 148: {"eitvbf", Pitvbf_e, -1},
! 149: {"eitvd", Pitvd_e, -1},
1.5 kondoh 150: #if 0
151: {"factorialitv",Pfactorialitv,1},
152: {"factorialitvd",Pfactorialitvd,1},
1.6 ! kondoh 153:
! 154: {"absitv", Pitv_abs, -2},
! 155: {"absitvbf", Pitvbf_abs, -2},
! 156: {"absitvd", Pitvd_abs, -2},
1.5 kondoh 157: #endif
158:
1.6 ! kondoh 159: {"logitv", Pitv_log, -2},
! 160: {"logitvbf", Pitvbf_log, -2},
! 161: {"logitvd", Pitvd_log, -2},
! 162: {"expitv", Pitv_exp, -2},
! 163: {"expitvbf", Pitvbf_exp, -2},
! 164: {"expitvd", Pitvd_exp, -2},
! 165: {"powitv", Pitv_pow, -3},
! 166: {"powitvbf", Pitvbf_pow, -3},
! 167: {"powitvd", Pitvd_pow, -3},
1.5 kondoh 168:
1.6 ! kondoh 169: {"sinitv", Pitv_sin, -2},
! 170: {"sinitvbf", Pitvbf_sin, -2},
1.5 kondoh 171: {"sinitvd", Pitvd_sin, -2},
1.6 ! kondoh 172: {"cositv", Pitv_cos, -2},
! 173: {"cositvbf", Pitvbf_cos, -2},
1.5 kondoh 174: {"cositvd", Pitvd_cos, -2},
1.6 ! kondoh 175: {"tanitv", Pitv_tan, -2},
! 176: {"tanitvbf", Pitvbf_tan, -2},
1.5 kondoh 177: {"tanitvd", Pitvd_tan, -2},
1.6 ! kondoh 178: {"asinitv", Pitv_asin, -2},
! 179: {"asinitvbf", Pitvbf_asin, -2},
! 180: {"asinitvd", Pitvd_asin, -2},
! 181: {"acositv", Pitv_acos, -2},
! 182: {"acositvbf", Pitvbf_acos, -2},
! 183: {"acositvd", Pitvd_acos, -2},
! 184: {"atanitv", Pitv_atan, -2},
! 185: {"atanitvbf", Pitvbf_atan, -2},
! 186: {"atanitvd", Pitvd_atan, -2},
! 187: {"sinhitv", Pitv_sinh, -2},
! 188: {"sinhitvbf", Pitvbf_sinh, -2},
! 189: {"sinhitvd", Pitvd_sinh, -2},
! 190: {"coshitv", Pitv_cosh, -2},
! 191: {"coshitvbf", Pitvbf_cosh, -2},
! 192: {"coshitvd", Pitvd_cosh, -2},
! 193: {"tanhitv", Pitv_tanh, -2},
! 194: {"tanhitvbf", Pitvbf_tanh, -2},
! 195: {"tanhitvd", Pitvd_tanh, -2},
! 196: {"asinhitv", Pitv_asinh, -2},
! 197: {"asinhitvbf", Pitvbf_asinh, -2},
! 198: {"asinhitvd", Pitvd_asinh, -2},
! 199: {"acoshitv", Pitv_acosh, -2},
! 200: {"acoshitvbf", Pitvbf_acosh, -2},
! 201: {"acoshitvd", Pitvd_acosh, -2},
! 202: {"atanhitv", Pitv_atanh, -2},
! 203: {"atanhitvbf", Pitvbf_atanh, -2},
! 204: {"atanhitvd", Pitvd_atanh, -2},
1.5 kondoh 205:
1.1 noro 206: /* plot time check */
207: {"ifcheck",Pifcheck,-7},
208: #endif
209: {0,0,0},
210: };
211:
1.3 kondoh 212: extern int mpfr_roundmode;
213:
1.1 noro 214: #if defined(INTERVAL)
215:
216: /* plot time check */
217: static void
218: Pifcheck(NODE arg, Obj *rp)
219: {
1.4 kondoh 220: Z m2,p2,s_id;
1.1 noro 221: NODE defrange;
222: LIST xrange,yrange,range[2],list,geom;
223: VL vl,vl0;
224: V v[2],av[2];
225: int ri,i,j,sign;
226: P poly;
227: P var;
228: NODE n,n0;
229: Obj t;
230:
231: struct canvas *can;
232: MAT m;
233: pointer **mb;
234: double **tabe, *px, *px1, *px2;
1.4 kondoh 235: Z one;
1.1 noro 236: int width, height, ix, iy;
237: int id;
238:
1.4 kondoh 239: STOZ(-2,m2); STOZ(2,p2);
240: STOZ(1,one);
1.1 noro 241: MKNODE(n,p2,0); MKNODE(defrange,m2,n);
242: poly = 0; vl = 0; geom = 0; ri = 0;
243: v[0] = v[1] = 0;
244: for ( ; arg; arg = NEXT(arg) ){
245: switch ( OID(BDY(arg)) ) {
246: case O_P:
247: poly = (P)BDY(arg);
248: get_vars_recursive((Obj)poly,&vl);
249: for(vl0=vl,i=0;vl0;vl0=NEXT(vl0)){
250: if(vl0->v->attr==(pointer)V_IND){
251: if(i>=2){
252: error("ifplot : invalid argument");
253: } else {
254: v[i++]=vl0->v;
255: }
256: }
257: }
258: break;
259: case O_LIST:
260: list = (LIST)BDY(arg);
261: if ( OID(BDY(BDY(list))) == O_P )
262: if ( ri > 1 )
263: error("ifplot : invalid argument");
264: else
265: range[ri++] = list;
266: else
267: geom = list;
268: break;
269: default:
270: error("ifplot : invalid argument"); break;
271: }
272: }
273: if ( !poly ) error("ifplot : invalid argument");
274: switch ( ri ) {
275: case 0:
276: if ( !v[1] ) error("ifplot : please specify all variables");
277: MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);
278: MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);
279: break;
280: case 1:
281: if ( !v[1] ) error("ifplot : please specify all variables");
282: av[0] = VR((P)BDY(BDY(range[0])));
283: if ( v[0] == av[0] ) {
284: xrange = range[0];
285: MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);
286: } else if ( v[1] == av[0] ) {
287: MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);
288: yrange = range[0];
289: } else
290: error("ifplot : invalid argument");
291: break;
292: case 2:
293: av[0] = VR((P)BDY(BDY(range[0])));
294: av[1] = VR((P)BDY(BDY(range[1])));
295: if ( ((v[0] == av[0]) && (!v[1] || v[1] == av[1])) ||
296: ((v[0] == av[1]) && (!v[1] || v[1] == av[0])) ) {
297: xrange = range[0]; yrange = range[1];
298: } else error("ifplot : invalid argument");
299: break;
300: default:
301: error("ifplot : cannot happen"); break;
302: }
303: can = canvas[id = search_canvas()];
304: if ( !geom ) {
305: width = 300;
306: height = 300;
307: can->width = 300;
308: can->height = 300;
309: } else {
1.4 kondoh 310: can->width = ZTOS((Z)BDY(BDY(geom)));
311: can->height = ZTOS((Z)BDY(NEXT(BDY(geom))));
1.1 noro 312: width = can->width;
313: height = can->height;
314: }
315: if ( xrange ) {
316: n = BDY(xrange); can->vx = VR((P)BDY(n)); n = NEXT(n);
317: can->qxmin = (Q)BDY(n); n = NEXT(n); can->qxmax = (Q)BDY(n);
318: can->xmin = ToReal(can->qxmin); can->xmax = ToReal(can->qxmax);
319: }
320: if ( yrange ) {
321: n = BDY(yrange); can->vy = VR((P)BDY(n)); n = NEXT(n);
322: can->qymin = (Q)BDY(n); n = NEXT(n); can->qymax = (Q)BDY(n);
323: can->ymin = ToReal(can->qymin); can->ymax = ToReal(can->qymax);
324: }
325: can->wname = "ifcheck";
326: can->formula = poly;
327: tabe = (double **)ALLOCA((width+1)*sizeof(double *));
328: for ( i = 0; i <= width; i++ )
329: tabe[i] = (double *)ALLOCA((height+1)*sizeof(double));
330: for(i=0;i<=width;i++)for(j=0;j<=height;j++)tabe[i][j]=0;
331: ccalc(tabe,can,0);
332: MKMAT(m,width,height);
333: mb = BDY(m);
334: for( ix=0; ix<width; ix++ ){
335: for( iy=0; iy<height; iy++){
336: if ( tabe[ix][iy] >= 0 ){
337: if ( (tabe[ix+1][iy] <= 0) ||
338: (tabe[ix][iy+1] <= 0 ) ||
339: (tabe[ix+1][iy+1] <= 0 ) ) mb[ix][iy] = (Obj)one;
340: } else {
341: if ( (tabe[ix+1][iy] >= 0 ) ||
342: ( tabe[ix][iy+1] >= 0 ) ||
343: ( tabe[ix+1][iy+1] >= 0 )) mb[ix][iy] = (Obj)one;
344: }
345: }
346: }
347: *rp = (Obj)m;
348: }
349:
350: void ccalc(double **tab,struct canvas *can,int nox)
351: {
352: double x,y,xmin,ymin,xstep,ystep;
353: int ix,iy;
354: Real r,rx,ry;
355: Obj fr,g;
356: int w,h;
357: V vx,vy;
358: Obj t,s;
359:
360: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
361: vx = can->vx;
362: vy = can->vy;
363: w = can->width; h = can->height;
364: xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
365: ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;
366: MKReal(1.0,rx); MKReal(1.0,ry);
367: for( ix = 0, x = xmin; ix < w+1 ; ix++, x += xstep ) {
368: BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
369: devalr(CO,t,&g);
370: for( iy = 0, y = ymin; iy < h+1 ; iy++, y += ystep ) {
371: BDY(ry) = y;
372: substr(CO,0,g,vy,y?(Obj)ry:0,&t);
373: devalr(CO,t,&s);
374: tab[ix][iy] = ToReal(s);
375: }
376: }
377: }
378: /* end plot time check */
379:
380: static void
1.3 kondoh 381: Pitvversion(NODE arg, Q *rp)
1.1 noro 382: {
1.4 kondoh 383: Z r;
384: STOZ(INT_ASIR_VERSION, r);
385: *rp = (Q)r;
1.1 noro 386: }
387:
388: extern int bigfloat;
389:
390: static void
391: Pitv(NODE arg, Obj *rp)
392: {
393: Num a, i, s;
394: Itv c;
395: double inf, sup;
396:
397: #if 1
398: if ( bigfloat )
399: Pitvbf(arg, rp);
400: else
401: Pitvd(arg,rp);
402: #else
403: asir_assert(ARG0(arg),O_N,"itv");
404: if ( argc(arg) > 1 ) {
405: asir_assert(ARG1(arg),O_N,"itv");
406: istoitv((Num)ARG0(arg),(Num)ARG1(arg),&c);
407: } else {
408: a = (Num)ARG0(arg);
409: if ( ! a ) {
410: *rp = 0;
411: return;
412: }
413: else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat) {
414: *rp = (Obj)a;
415: return;
416: }
417: else if ( NID(a) == N_IntervalDouble ) {
418: inf = INF((IntervalDouble)a);
419: sup = SUP((IntervalDouble)a);
420: double2bf(inf, (BF *)&i);
421: double2bf(sup, (BF *)&s);
422: istoitv(i,s,&c);
423: }
424: else istoitv(a,a,&c);
425: }
426: if ( NID( c ) == N_IntervalBigFloat )
427: addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);
428: else *rp = (Obj)c;
429: #endif
430: }
431:
432: static void
433: Pitvbf(NODE arg, Obj *rp)
434: {
435: Num a, i, s;
1.3 kondoh 436: IntervalBigFloat c;
437: Num ii,ss;
438: Real di, ds;
1.1 noro 439: double inf, sup;
1.3 kondoh 440: int current_roundmode;
1.1 noro 441:
442: asir_assert(ARG0(arg),O_N,"intvalbf");
443: a = (Num)ARG0(arg);
444: if ( argc(arg) > 1 ) {
445: asir_assert(ARG1(arg),O_N,"intvalbf");
1.3 kondoh 446:
1.1 noro 447: i = (Num)ARG0(arg);
448: s = (Num)ARG1(arg);
1.3 kondoh 449: current_roundmode = mpfr_roundmode;
450: mpfr_roundmode = MPFR_RNDD;
451: ii = tobf(i, DEFAULTPREC);
452: mpfr_roundmode = MPFR_RNDU;
453: ss = tobf(s, DEFAULTPREC);
454: istoitv(ii,ss,(Itv *)&c);
455: // MKIntervalBigFloat((BF)ii,(BF)ss,c);
456: // ToBf(s, &ss);
457: mpfr_roundmode = current_roundmode;
1.1 noro 458: } else {
459: if ( ! a ) {
460: *rp = 0;
461: return;
462: }
463: else if ( NID(a) == N_IP ) {
464: itvtois((Itv)a, &i, &s);
1.3 kondoh 465: current_roundmode = mpfr_roundmode;
466: mpfr_roundmode = MPFR_RNDD;
467: ii = tobf(i, DEFAULTPREC);
468: mpfr_roundmode = MPFR_RNDU;
469: ss = tobf(s, DEFAULTPREC);
470: istoitv(ii,ss,(Itv *)&c);
471: // MKIntervalBigFloat((BF)ii,(BF)ss,c);
472: mpfr_roundmode = current_roundmode;
1.1 noro 473: }
474: else if ( NID(a) == N_IntervalBigFloat) {
475: *rp = (Obj)a;
476: return;
477: }
478: else if ( NID(a) == N_IntervalDouble ) {
479: inf = INF((IntervalDouble)a);
480: sup = SUP((IntervalDouble)a);
1.3 kondoh 481: current_roundmode = mpfr_roundmode;
482: //double2bf(inf, (BF *)&i);
483: //double2bf(sup, (BF *)&s);
484: mpfr_roundmode = MPFR_RNDD;
485: MKReal(inf,di);
486: ii = tobf((Num)di, DEFAULTPREC);
487: mpfr_roundmode = MPFR_RNDU;
488: MKReal(sup,ds);
489: ss = tobf((Num)ds, DEFAULTPREC);
490: istoitv(ii,ss,(Itv *)&c);
491: // MKIntervalBigFloat((BF)ii,(BF)ss,c);
492: mpfr_roundmode = current_roundmode;
1.1 noro 493: }
494: else {
1.3 kondoh 495: current_roundmode = mpfr_roundmode;
496: mpfr_roundmode = MPFR_RNDD;
497: ii = tobf(a, DEFAULTPREC);
498: mpfr_roundmode = MPFR_RNDU;
499: ss = tobf(a, DEFAULTPREC);
500: //ToBf(a, (BF *)&i);
501: istoitv(ii,ss,(Itv *)&c);
502: // MKIntervalBigFloat((BF)ii,(BF)ss,c);
503: mpfr_roundmode = current_roundmode;
1.1 noro 504: }
505: }
1.3 kondoh 506: // if ( c && OID( c ) == O_N && NID( c ) == N_IntervalBigFloat )
507: // addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);
508: // else *rp = (Obj)c;
509: *rp = (Obj)c;
1.1 noro 510: }
511:
512: static void
513: Pitvd(NODE arg, Obj *rp)
514: {
515: double inf, sup;
516: Num a, a0, a1, t;
517: Itv ia;
518: IntervalDouble d;
519:
520: asir_assert(ARG0(arg),O_N,"intvald");
521: a0 = (Num)ARG0(arg);
522: if ( argc(arg) > 1 ) {
523: asir_assert(ARG1(arg),O_N,"intvald");
524: a1 = (Num)ARG1(arg);
525: } else {
526: if ( a0 && OID(a0)==O_N && NID(a0)==N_IntervalDouble ) {
527: inf = INF((IntervalDouble)a0);
528: sup = SUP((IntervalDouble)a0);
529: MKIntervalDouble(inf,sup,d);
530: *rp = (Obj)d;
531: return;
532: }
533: a1 = (Num)ARG0(arg);
534: }
535: if ( compnum(0,a0,a1) > 0 ) {
536: t = a0; a0 = a1; a1 = t;
537: }
1.5 kondoh 538: inf = toRealDown(a0);
539: sup = toRealUp(a1);
1.1 noro 540: MKIntervalDouble(inf,sup,d);
541: *rp = (Obj)d;
542: }
543:
544: static void
545: Pinf(NODE arg, Obj *rp)
546: {
547: Num a, i, s;
548: Real r;
549: double d;
550:
551: a = (Num)ARG0(arg);
552: if ( ! a ) {
553: *rp = 0;
554: } else if ( OID(a) == O_N ) {
555: switch ( NID(a) ) {
556: case N_IntervalDouble:
557: d = INF((IntervalDouble)a);
558: MKReal(d, r);
559: *rp = (Obj)r;
560: break;
561: case N_IP:
562: case N_IntervalBigFloat:
563: case N_IntervalQuad:
564: itvtois((Itv)ARG0(arg),&i,&s);
565: *rp = (Obj)i;
566: break;
567: default:
568: *rp = (Obj)a;
569: break;
570: }
571: } else {
572: *rp = (Obj)a;
573: }
574: }
575:
576: static void
577: Psup(NODE arg, Obj *rp)
578: {
579: Num a, i, s;
580: Real r;
581: double d;
582:
583: a = (Num)ARG0(arg);
584: if ( ! a ) {
585: *rp = 0;
586: } else if ( OID(a) == O_N ) {
587: switch ( NID(a) ) {
588: case N_IntervalDouble:
589: d = SUP((IntervalDouble)a);
590: MKReal(d, r);
591: *rp = (Obj)r;
592: break;
593: case N_IP:
594: case N_IntervalBigFloat:
595: case N_IntervalQuad:
596: itvtois((Itv)ARG0(arg),&i,&s);
597: *rp = (Obj)s;
598: break;
599: default:
600: *rp = (Obj)a;
601: break;
602: }
603: } else {
604: *rp = (Obj)a;
605: }
606: }
607:
608: static void
609: Pmid(NODE arg, Obj *rp)
610: {
611: Num a, s;
612: Real r;
613: double d;
614:
615: a = (Num)ARG0(arg);
616: if ( ! a ) {
617: *rp = 0;
618: } else switch (OID(a)) {
619: case O_N:
620: if ( NID(a) == N_IntervalDouble ) {
621: d = ( INF((IntervalDouble)a)+SUP((IntervalDouble)a) ) / 2.0;
622: MKReal(d, r);
623: *rp = (Obj)r;
624: } else if ( NID(a) == N_IntervalQuad ) {
625: error("mid: not supported operation");
626: *rp = 0;
627: } else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat ) {
628: miditvp((Itv)ARG0(arg),&s);
629: *rp = (Obj)s;
630: } else {
631: *rp = (Obj)a;
632: }
633: break;
634: #if 0
635: case O_P:
636: case O_R:
637: case O_LIST:
638: case O_VECT:
639: case O_MAT:
640: #endif
641: default:
642: *rp = (Obj)a;
643: break;
644: }
645: }
646:
647: static void
648: Pcup(NODE arg, Obj *rp)
649: {
650: Itv s;
651: Num a, b;
652:
653: asir_assert(ARG0(arg),O_N,"cup");
654: asir_assert(ARG1(arg),O_N,"cup");
655: a = (Num)ARG0(arg);
656: b = (Num)ARG1(arg);
657: if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
658: cupitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);
659: } else {
660: cupitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
661: *rp = (Obj)s;
662: }
663: }
664:
665: static void
666: Pcap(NODE arg, Obj *rp)
667: {
668: Itv s;
669: Num a, b;
670:
671: asir_assert(ARG0(arg),O_N,"cap");
672: asir_assert(ARG1(arg),O_N,"cap");
673: a = (Num)ARG0(arg);
674: b = (Num)ARG1(arg);
675: if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
676: capitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);
677: } else {
678: capitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
679: *rp = (Obj)s;
680: }
681: }
682:
683: static void
684: Pwidth(arg,rp)
685: NODE arg;
686: Obj *rp;
687: {
688: Num s;
689: Num a;
690:
691: asir_assert(ARG0(arg),O_N,"width");
692: a = (Num)ARG0(arg);
693: if ( ! a ) {
694: *rp = 0;
695: } else if ( NID(a) == N_IntervalDouble ) {
696: widthitvd((IntervalDouble)a, (Num *)rp);
697: } else {
698: widthitvp((Itv)ARG0(arg),&s);
699: *rp = (Obj)s;
700: }
701: }
702:
703: static void
704: Pabsitv(arg,rp)
705: NODE arg;
706: Obj *rp;
707: {
708: Num s;
709: Num a, b;
710:
711: asir_assert(ARG0(arg),O_N,"absitv");
712: a = (Num)ARG0(arg);
713: if ( ! a ) {
714: *rp = 0;
715: } else if ( NID(a) == N_IntervalDouble ) {
716: absitvd((IntervalDouble)a, (Num *)rp);
717: } else {
718: absitvp((Itv)ARG0(arg),&s);
719: *rp = (Obj)s;
720: }
721: }
722:
723: static void
724: Pdistance(arg,rp)
725: NODE arg;
726: Obj *rp;
727: {
728: Num s;
729: Num a, b;
730:
731: asir_assert(ARG0(arg),O_N,"distance");
732: asir_assert(ARG1(arg),O_N,"distance");
733: a = (Num)ARG0(arg);
734: b = (Num)ARG1(arg);
735: if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
736: distanceitvd((IntervalDouble)a, (IntervalDouble)b, (Num *)rp);
737: } else {
738: distanceitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
739: *rp = (Obj)s;
740: }
741: }
742:
743: static void
1.4 kondoh 744: Pinitv(NODE arg, Obj *rp)
1.1 noro 745: {
746: int s;
1.4 kondoh 747: Z q;
1.1 noro 748:
749: asir_assert(ARG0(arg),O_N,"intval");
750: asir_assert(ARG1(arg),O_N,"intval");
751: if ( ! ARG1(arg) ) {
752: if ( ! ARG0(arg) ) s = 1;
753: else s = 0;
754: }
755: else if ( NID(ARG1(arg)) == N_IntervalDouble ) {
756: s = initvd((Num)ARG0(arg),(IntervalDouble)ARG1(arg));
757:
758: } else if ( NID(ARG1(arg)) == N_IP || NID(ARG1(arg)) == N_IntervalBigFloat ) {
759: if ( ! ARG0(arg) ) s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));
760: else if ( NID(ARG0(arg)) == N_IP ) {
761: s = itvinitvp((Itv)ARG0(arg),(Itv)ARG1(arg));
762: } else {
763: s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));
764: }
765: } else {
766: s = ! compnum(0,(Num)ARG0(arg),(Num)ARG1(arg));
767: }
1.4 kondoh 768: STOZ(s,q);
1.1 noro 769: *rp = (Obj)q;
770: }
771:
772: static void
773: Pdisjitv(arg,rp)
774: NODE arg;
775: Obj *rp;
776: {
777: Itv s;
778:
779: asir_assert(ARG0(arg),O_N,"disjitv");
780: asir_assert(ARG1(arg),O_N,"disjitv");
781: error("disjitv: not implemented yet");
782: if ( ! s ) *rp = 0;
783: else *rp = (Obj)ONE;
784: }
785:
1.3 kondoh 786: static void
787: PzeroRewriteMode(NODE arg, Obj *rp)
788: {
1.4 kondoh 789: Q a;
790: Z r;
1.3 kondoh 791:
1.4 kondoh 792: STOZ(zerorewrite,r);
1.3 kondoh 793: *rp = (Obj)r;
794:
795: if (arg) {
796: a = (Q)ARG0(arg);
797: if(!a) {
798: zerorewrite = 0;
799: } else if ( (NUM(a)&&INT(a)) ){
800: zerorewrite = 1;
801: }
802: }
803: }
804:
805: static void
806: PzeroRewriteCountClear(NODE arg, Obj *rp)
807: {
1.4 kondoh 808: Q a;
809: Z r;
1.3 kondoh 810:
1.4 kondoh 811: STOZ(zerorewriteCount,r);
1.3 kondoh 812: *rp = (Obj)r;
813:
814: if (arg) {
815: a = (Q)ARG0(arg);
816: if(a &&(NUM(a)&&INT(a))){
817: zerorewriteCount = 0;
818: }
819: }
820: }
821:
822: static void
823: PzeroRewriteCount(NODE arg, Obj *rp)
824: {
1.4 kondoh 825: Z r;
1.3 kondoh 826:
1.4 kondoh 827: STOZ(zerorewriteCount,r);
1.3 kondoh 828: *rp = (Obj)r;
829: }
830:
831:
1.1 noro 832: #endif
833: extern int printmode;
834:
835: static void pprintmode( void )
836: {
837: switch (printmode) {
838: #if defined(INTERVAL)
839: case MID_PRINTF_E:
840: fprintf(stderr,"Interval printing mode is a mitpoint type.\n");
841: #endif
842: case PRINTF_E:
843: fprintf(stderr,"Printf's double printing mode is \"%%.16e\".\n");
844: break;
845: #if defined(INTERVAL)
846: case MID_PRINTF_G:
847: fprintf(stderr,"Interval printing mode is a mitpoint type.\n");
848: #endif
849: default:
850: case PRINTF_G:
851: fprintf(stderr,"Printf's double printing mode is \"%%g\".\n");
852: break;
853: }
854: }
855:
856: static void
857: Pprintmode(NODE arg, Obj *rp)
858: {
859: int l;
1.4 kondoh 860: Z a, r;
1.1 noro 861:
1.4 kondoh 862: a = (Z)ARG0(arg);
1.1 noro 863: if(!a||(NUM(a)&&INT(a))){
1.4 kondoh 864: l=ZTOS(a);
1.1 noro 865: if ( l < 0 ) l = 0;
866: #if defined(INTERVAL)
867: else if ( l > MID_PRINTF_E ) l = 0;
868: #else
869: else if ( l > PRINTF_E ) l = 0;
870: #endif
1.4 kondoh 871: STOZ(printmode,r);
1.1 noro 872: *rp = (Obj)r;
873: printmode = l;
874: pprintmode();
875: } else {
876: *rp = 0;
877: }
878: }
1.5 kondoh 879:
880: #if defined(INTERVAL)
881: void
882: Ppi_itvd(NODE arg, Obj *rp)
883: {
884: double inf, sup;
885: IntervalDouble c;
886: FPMINUSINF
887: sscanf("3.1415926535897932384626433832795028841971693993751", "%lf", &inf);
888: FPPLUSINF
889: sscanf("3.1415926535897932384626433832795028841971693993752", "%lf", &sup);
890: FPNEAREST
891: MKIntervalDouble(inf,sup,c);
892: *rp = (Obj)c;
893: }
894: void
895: Pe_itvd(NODE arg, Obj *rp)
896: {
897: double inf, sup;
898: IntervalDouble c;
899: FPMINUSINF
900: sscanf( "2.7182818284590452353602874713526624977572470936999", "%lf", &inf);
901: FPPLUSINF
902: sscanf( "2.7182818284590452353602874713526624977572470937000", "%lf", &sup);
903: FPNEAREST
904: MKIntervalDouble(inf,sup,c);
905: *rp = (Obj)c;
906: }
907: void
908: Pln2_itvd(NODE arg, Obj *rp)
909: {
910: double inf, sup;
911: IntervalDouble c;
912: FPMINUSINF
913: sscanf( "0.69314718055994530941723212145817656807550013436025", "%lf", &inf);
914: FPPLUSINF
915: sscanf( "0.69314718055994530941723212145817656807550013436026", "%lf", &sup);
916: FPNEAREST
917: MKIntervalDouble(inf,sup,c);
918: *rp = (Obj)c;
919: }
920:
921: void mpfi_func(NODE arg, int (*mpfi_f)(), int prec, Obj *rp)
922: {
923: Num a, ii, ss;
924: Itv c;
925: BF inf, sup;
926: int arg1prec;
927: mpfi_t mpitv, rv;
928:
929:
930: /*
931: if ( argc(arg) == 2 ) {
932: prec = QTOS((Q)ARG1(arg))*3.32193;
933: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
934: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
935: } else {
936: prec = 0;
937: prec = mpfr_get_default_prec();
938: }
939: */
940: if ( prec > 0 ) arg1prec = prec;
941: else arg1prec = NEXT(arg) ? ZTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
942: a = ARG0(arg);
943: itvtois((Itv)a, &ii, &ss);
944:
945: inf = (BF)tobf(ii, arg1prec);
946: sup = (BF)tobf(ss, arg1prec);
947:
948: mpfi_init2(rv,arg1prec);
949: mpfi_init2(mpitv,arg1prec);
950: mpfr_set(&(mpitv->left), BDY(inf), MPFR_RNDD);
951: mpfr_set(&(mpitv->right), BDY(sup), MPFR_RNDU);
952:
953: (*mpfi_f)(rv, mpitv);
954: //mpfi_sin(rv, mpitv);
955:
956: MPFRTOBF(&(rv->left), inf);
957: MPFRTOBF(&(rv->right), sup);
958:
959: if ( !cmpbf((Num)inf,0) ) inf = 0;
960: if ( !cmpbf((Num)sup,0) ) sup = 0;
961:
962: if ( inf || sup ) {
963: istoitv((Num)inf, (Num)sup, &c);
964: } else {
965: c = 0;
966: }
967: *rp = (Obj)c;
968: //mpfi_clear(rv);
969: mpfi_clear(mpitv);
970: }
971:
972: void mpfi_func_d(NODE arg, int (*mpfi_f)(), Obj *rp)
973: {
974: Obj bfv;
975: Num ii, ss;
976: IntervalDouble c;
977: double inf, sup;
978:
979: mpfi_func(arg, mpfi_f, 53, &bfv);
980: itvtois((Itv)bfv, &ii, &ss);
981: inf = toRealDown(ii);
982: sup = toRealUp(ss);
983: MKIntervalDouble(inf,sup,c);
984: *rp = (Obj)c;
985: }
986:
1.1 noro 987:
1.5 kondoh 988: void
989: Psinitvd(NODE arg, Obj *rp)
990: {
991: Obj bfv;
992: Num ii,ss;
993: IntervalDouble c;
994: double ai,as,mas, bi,bs;
995: double inf,sup;
996:
997: mpfi_func(arg, mpfi_sin, 53, &bfv);
998: itvtois((Itv)bfv, &ii, &ss);
999: inf = toRealDown(ii);
1000: sup = toRealUp(ss);
1001: MKIntervalDouble(inf,sup,c);
1002: *rp = (Obj)c;
1003: }
1004:
1.6 ! kondoh 1005: static
1.5 kondoh 1006: void
1007: Psinitv(NODE arg, Obj *rp)
1008: {
1009: mpfi_func(arg, mpfi_sin, 0, rp);
1010: }
1011:
1012:
1013:
1014:
1015: //void evalitvr(VL, Obj, int, int, Obj *);
1016:
1017: static void
1018: Pevalitv(NODE arg, Obj *rp)
1019: {
1.6 ! kondoh 1020: if ( bigfloat )
! 1021: Pevalitvbf(arg, rp);
! 1022: else
! 1023: Pevalitvd(arg, rp);
! 1024: }
! 1025:
! 1026: static void
! 1027: Pevalitvbf(NODE arg, Obj *rp)
! 1028: {
1.5 kondoh 1029: int prec;
1030:
1031: asir_assert(ARG0(arg),O_R,"evalitv");
1032: if ( argc(arg) == 2 ) {
1033: long int mpfr_prec_max = MPFR_PREC_MAX;
1034: prec = ZTOS((Q)ARG1(arg))*3.32193;
1035: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
1036: //else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
1037: else if ( prec > mpfr_prec_max ) prec = mpfr_prec_max;
1038: } else
1039: prec = 0;
1040: evalitvr(CO,(Obj)ARG0(arg),prec,EvalIntervalBigFloat,rp);
1041: }
1042:
1043: static void
1044: Pevalitvd(NODE arg, Obj *rp)
1045: {
1046: asir_assert(ARG0(arg),O_R,"evalitvd");
1047: evalitvr(CO,(Obj)ARG0(arg),53,EvalIntervalDouble,rp);
1048: }
1049:
1050: // in parse/puref.c
1051: void instoobj(PFINS ins,Obj *rp);
1052:
1053: // in this
1054: void evalitvr(VL, Obj, int, int, Obj *);
1055: void evalitvp(VL, P, int, int, P *);
1056: void evalitvv(VL, V, int, int, Obj *);
1057: void evalitvins(PFINS, int, int, Obj *);
1058:
1059:
1060:
1061: void evalitvr(VL vl,Obj a,int prec, int type, Obj *c)
1062: {
1063: Obj nm,dn;
1064:
1065: if ( !a )
1066: *c = 0;
1067: else {
1068: switch ( OID(a) ) {
1069: case O_N:
1070: toInterval((Num)a, prec, type, (Num *)c);
1071: break;
1072: case O_P:
1073: evalitvp(vl,(P)a,prec,type,(P *)c);
1074: break;
1075: case O_R:
1076: evalitvp(vl,NM((R)a),prec,type,(P *)&nm);
1077: evalitvp(vl,DN((R)a),prec,type,(P *)&dn);
1078: divr(vl,nm,dn,c);
1079: break;
1080: default:
1081: error("evalr : not implemented"); break;
1082: }
1083: }
1084: }
1085:
1086: void evalitvp(VL vl,P p,int prec, int type, P *pr)
1087: {
1088: P t;
1089: DCP dc,dcr0,dcr;
1090: Obj u;
1091:
1092: if ( !p || NUM(p) ) {
1093: toInterval((Num)p, prec, type, (Num *)pr);
1094: //*pr = p;
1095: } else {
1096: for ( dcr0 = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
1097: evalitvp(vl,COEF(dc),prec,type, &t);
1098: if ( t ) {
1099: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
1100: }
1101: }
1102: if ( !dcr0 ) {
1103: *pr = 0; return;
1104: } else {
1105: NEXT(dcr) = 0; MKP(VR(p),dcr0,t);
1106: }
1107: if ( NUM(t) ) {
1108: //*pr = t;
1109: toInterval((Num)t, prec, type, (Num *)pr);
1110: return;
1111: } else if ( (VR(t) != VR(p)) || ((vid)VR(p)->attr != V_PF) ) {
1112: *pr = t; return;
1113: } else {
1114: evalitvv(vl,VR(p),prec,type,&u);
1115: substr(vl,1,(Obj)t,VR(p),u, (Obj *)&t);
1116: if ( t && NUM(t) ) {
1117: toInterval((Num)t, prec, type, (Num *)pr);
1118: } else {
1119: *pr = t;
1120: }
1121: }
1122: }
1123: }
1124:
1125: void evalitvv(VL vl,V v,int prec, int type, Obj *rp)
1126: {
1127: PFINS ins,tins;
1128: PFAD ad,tad;
1129: PF pf;
1130: P t;
1131: int i;
1132:
1133: if ( (vid)v->attr != V_PF ) {
1134: MKV(v,t); *rp = (Obj)t;
1135: } else {
1136: ins = (PFINS)v->priv; ad = ins->ad; pf = ins->pf;
1137: tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
1138: tins->pf = pf;
1139: for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) {
1140: tad[i].d = ad[i].d;
1141: evalitvr(vl,ad[i].arg,prec,type,&tad[i].arg);
1142: }
1143: evalitvins(tins,prec,type,rp);
1144: }
1145: }
1146:
1147: void evalitvins(PFINS ins,int prec, int type, Obj *rp)
1148: {
1149: PF pf;
1150: PFINS tins;
1151: PFAD ad,tad;
1152: int i;
1153: Z q;
1154: V v;
1155: P x;
1156: NODE n0,n;
1157:
1158:
1159: pf = ins->pf; ad = ins->ad;
1160: tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
1161: tins->pf = pf; tad = tins->ad;
1162: for ( i = 0; i < pf->argc; i++ ) {
1163: //tad[i].d = ad[i].d; evalr(CO,ad[i].arg,prec,&tad[i].arg);
1164: tad[i].d = ad[i].d; evalitvr(CO,ad[i].arg,prec,type,&tad[i].arg);
1165: }
1166: for ( i = 0; i < pf->argc; i++ )
1167: if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break;
1168: if ( (i != pf->argc) || !pf->intervalfunc[type] ) { ///////////////////////////
1169: instoobj(tins,rp);
1170: } else {
1171: int IsCPLX = 0;
1172: for ( n0 = 0, i = 0; i < pf->argc; i++ ) {
1173: NEXTNODE(n0,n); BDY(n) = (pointer)tad[i].arg;
1174: if (tad[i].arg && NID(tad[i].arg) == N_C) IsCPLX = 1;
1175: }
1176: if ( prec ) {
1177: NEXTNODE(n0,n); STOZ(prec,q); BDY(n) = (pointer)q;
1178: }
1179: if ( n0 ) NEXT(n) = 0;
1180:
1181:
1182: if ( IsCPLX ) {
1183: instoobj(tins,rp);
1184: } else {
1185: (*pf->intervalfunc[type])(n0,rp);
1186: }
1187: }
1188: }
1189:
1190:
1.6 ! kondoh 1191: static
1.5 kondoh 1192: void Pitvbf_pi(NODE arg, Obj *rp)
1193: {
1194: BF inf, sup;
1195: IntervalBigFloat c;
1196: mpfi_t rv;
1197: int prec;
1198:
1199: prec = (arg) ? ZTOS((Q)ARG0(arg)) : 0; //mpfr_get_default_prec();
1200: prec ? mpfi_init2(rv,prec) : mpfi_init(rv);
1201:
1202: mpfi_const_pi(rv);
1203:
1204: MPFRTOBF(&(rv->left), inf);
1205: MPFRTOBF(&(rv->right), sup);
1206:
1207: MKIntervalBigFloat(inf,sup,c);
1208:
1209: *rp = (Obj)c;
1210: }
1.1 noro 1211:
1.6 ! kondoh 1212: static
1.5 kondoh 1213: void Pitvd_pi(NODE arg, Obj *rp)
1214: {
1215: BF bfinf, bfsup;
1216: IntervalDouble c;
1217: mpfi_t rv;
1218: double inf, sup;
1219:
1220: mpfi_init2(rv,53);
1221:
1222: mpfi_const_pi(rv);
1223:
1224: MPFRTOBF(&(rv->left), bfinf);
1225: MPFRTOBF(&(rv->right), bfsup);
1226:
1227: inf = toRealDown((Num)bfinf);
1228: sup = toRealUp((Num)bfsup);
1229:
1230: MKIntervalDouble(inf,sup,c);
1231:
1232: *rp = (Obj)c;
1233: }
1234:
1.6 ! kondoh 1235: static
1.5 kondoh 1236: void Pitvbf_e(NODE arg,Obj *rp)
1237: {
1238: BF inf, sup;
1239: IntervalBigFloat c;
1240: mpfi_t rv;
1241: mpfi_t one;
1242: int prec;
1243:
1244: prec = (arg) ? ZTOS((Q)ARG0(arg)) : 0; //mpfr_get_default_prec();
1245: prec ? mpfi_init2(rv,prec) : mpfi_init(rv);
1246:
1247: mpfi_init(one);
1248: mpfi_set_ui(one,1);
1249: mpfi_exp(rv,one);
1250:
1251: MPFRTOBF(&(rv->left), inf);
1252: MPFRTOBF(&(rv->right), sup);
1253:
1254: MKIntervalBigFloat(inf,sup,c);
1255: //istoitv((Num)inf, (Num)sup, &c);
1256: *rp = (Obj)c;
1257: mpfi_clear(one);
1258: }
1259:
1.6 ! kondoh 1260: static
1.5 kondoh 1261: void Pitvd_e(NODE arg, Obj *rp)
1262: {
1263: BF bfinf, bfsup;
1264: IntervalDouble c;
1265: mpfi_t rv;
1266: mpfi_t one;
1267: double inf, sup;
1268:
1269: mpfi_init2(rv,53);
1270:
1271: mpfi_init2(one, 53);
1272: mpfi_set_ui(one,1);
1273: mpfi_exp(rv,one);
1274:
1275:
1276: MPFRTOBF(&(rv->left), bfinf);
1277: MPFRTOBF(&(rv->right), bfsup);
1278:
1279: inf = toRealDown((Num)bfinf);
1280: sup = toRealUp((Num)bfsup);
1281:
1282: MKIntervalDouble(inf,sup,c);
1283:
1284: *rp = (Obj)c;
1285: }
1286:
1287: void (*pi_itv_ft[])() = {Pitvd_pi, 0, Pitvbf_pi};
1288: void (*e_itv_ft[])() = {Pitvd_e, 0, Pitvbf_e};
1289: //void (*sin_itv_ft[])() = {Psinitvd, 0, Psinitv};
1290: void (*sin_itv_ft[])() = {Pitvd_sin, 0, Pitvbf_sin};
1291: void (*cos_itv_ft[])() = {Pitvbf_cos, 0, Pitvbf_cos};
1292: void (*tan_itv_ft[])() = {Pitvd_tan, 0, Pitvbf_tan};
1293: void (*asin_itv_ft[])() = {Pitvd_asin, 0, Pitvbf_asin};
1294: void (*acos_itv_ft[])() = {Pitvd_acos, 0, Pitvbf_acos};
1295: void (*atan_itv_ft[])() = {Pitvd_atan, 0, Pitvbf_atan};
1296: void (*sinh_itv_ft[])() = {Pitvd_sinh, 0, Pitvbf_sinh};
1297: void (*cosh_itv_ft[])() = {Pitvd_cosh, 0, Pitvbf_cosh};
1298: void (*tanh_itv_ft[])() = {Pitvd_tanh, 0, Pitvbf_tanh};
1299: void (*asinh_itv_ft[])() = {Pitvd_asinh, 0, Pitvbf_asinh};
1300: void (*acosh_itv_ft[])() = {Pitvd_acosh, 0, Pitvbf_acosh};
1301: void (*atanh_itv_ft[])() = {Pitvd_atanh, 0, Pitvbf_atanh};
1302: void (*exp_itv_ft[])() = {Pitvd_exp, 0, Pitvbf_exp};
1303: void (*log_itv_ft[])() = {Pitvd_log, 0, Pitvbf_log};
1304: void (*abs_itv_ft[])() = {0};
1305: void (*pow_itv_ft[])() = {Pitvbf_pow, 0, Pitvbf_pow};
1306: //void (*pow_itv_ft[])() = {0, 0, 0};
1307:
1308:
1.6 ! kondoh 1309: static
1.5 kondoh 1310: void Pitvbf_sin(NODE arg,Obj *rp)
1311: {
1312: mpfi_func(arg, mpfi_sin, 0, rp);
1313: }
1314:
1.6 ! kondoh 1315: static
1.5 kondoh 1316: void Pitvbf_cos(NODE arg,Obj *rp)
1317: {
1318: mpfi_func(arg, mpfi_cos, 0, rp);
1319: }
1320:
1.6 ! kondoh 1321: static
1.5 kondoh 1322: void Pitvbf_tan(NODE arg,Obj *rp)
1323: {
1324: mpfi_func(arg, mpfi_tan, 0, rp);
1325: }
1326:
1.6 ! kondoh 1327: static
1.5 kondoh 1328: void Pitvbf_asin(NODE arg,Obj *rp)
1329: {
1330: mpfi_func(arg, mpfi_asin, 0, rp);
1331: }
1332:
1.6 ! kondoh 1333: static
1.5 kondoh 1334: void Pitvbf_acos(NODE arg,Obj *rp)
1335: {
1336: mpfi_func(arg, mpfi_acos, 0, rp);
1337: }
1338:
1.6 ! kondoh 1339: static
1.5 kondoh 1340: void Pitvbf_atan(NODE arg,Obj *rp)
1341: {
1342: mpfi_func(arg, mpfi_atan, 0, rp);
1343: }
1344:
1.6 ! kondoh 1345: static
1.5 kondoh 1346: void Pitvbf_sinh(NODE arg,Obj *rp)
1347: {
1348: mpfi_func(arg, mpfi_sinh, 0, rp);
1349: }
1350:
1.6 ! kondoh 1351: static
1.5 kondoh 1352: void Pitvbf_cosh(NODE arg,Obj *rp)
1353: {
1354: mpfi_func(arg, mpfi_cosh, 0, rp);
1355: }
1356:
1.6 ! kondoh 1357: static
1.5 kondoh 1358: void Pitvbf_tanh(NODE arg,Obj *rp)
1359: {
1360: mpfi_func(arg, mpfi_tanh, 0, rp);
1361: }
1362:
1.6 ! kondoh 1363: static
1.5 kondoh 1364: void Pitvbf_asinh(NODE arg,Obj *rp)
1365: {
1366: mpfi_func(arg, mpfi_asinh, 0, rp);
1367: }
1368:
1.6 ! kondoh 1369: static
1.5 kondoh 1370: void Pitvbf_acosh(NODE arg,Obj *rp)
1371: {
1372: mpfi_func(arg, mpfi_acosh, 0, rp);
1373: }
1374:
1.6 ! kondoh 1375: static
1.5 kondoh 1376: void Pitvbf_atanh(NODE arg,Obj *rp)
1377: {
1378: mpfi_func(arg, mpfi_atanh, 0, rp);
1379: }
1380:
1.6 ! kondoh 1381: static
1.5 kondoh 1382: void Pitvbf_exp(NODE arg,Obj *rp)
1383: {
1384: mpfi_func(arg, mpfi_exp, 0, rp);
1385: }
1386:
1.6 ! kondoh 1387: static
1.5 kondoh 1388: void Pitvbf_log(NODE arg,Obj *rp)
1389: {
1390: mpfi_func(arg, mpfi_log, 0, rp);
1391: }
1392:
1.6 ! kondoh 1393: static
1.5 kondoh 1394: void Pitvbf_abs(NODE arg,Obj *rp)
1395: {
1396: mpfi_func(arg, mpfi_abs, 0, rp);
1397: }
1398:
1.6 ! kondoh 1399: static
1.5 kondoh 1400: void Pitvd_sin(NODE arg,Obj *rp)
1401: {
1402: mpfi_func_d(arg, mpfi_sin, rp);
1403: }
1404:
1.6 ! kondoh 1405: static
1.5 kondoh 1406: void Pitvd_cos(NODE arg,Obj *rp)
1407: {
1408: mpfi_func_d(arg, mpfi_cos, rp);
1409: }
1410:
1.6 ! kondoh 1411: static
1.5 kondoh 1412: void Pitvd_tan(NODE arg,Obj *rp)
1413: {
1414: mpfi_func_d(arg, mpfi_tan, rp);
1415: }
1416:
1.6 ! kondoh 1417: static
1.5 kondoh 1418: void Pitvd_asin(NODE arg,Obj *rp)
1419: {
1420: mpfi_func_d(arg, mpfi_asin, rp);
1421: }
1422:
1.6 ! kondoh 1423: static
1.5 kondoh 1424: void Pitvd_acos(NODE arg,Obj *rp)
1425: {
1426: mpfi_func_d(arg, mpfi_acos, rp);
1427: }
1428:
1.6 ! kondoh 1429: static
1.5 kondoh 1430: void Pitvd_atan(NODE arg,Obj *rp)
1431: {
1432: mpfi_func_d(arg, mpfi_atan, rp);
1433: }
1434:
1.6 ! kondoh 1435: static
1.5 kondoh 1436: void Pitvd_sinh(NODE arg,Obj *rp)
1437: {
1438: mpfi_func_d(arg, mpfi_sinh, rp);
1439: }
1440:
1.6 ! kondoh 1441: static
1.5 kondoh 1442: void Pitvd_cosh(NODE arg,Obj *rp)
1443: {
1444: mpfi_func_d(arg, mpfi_cosh, rp);
1445: }
1446:
1.6 ! kondoh 1447: static
1.5 kondoh 1448: void Pitvd_tanh(NODE arg,Obj *rp)
1449: {
1450: mpfi_func_d(arg, mpfi_tanh, rp);
1451: }
1452:
1.6 ! kondoh 1453: static
1.5 kondoh 1454: void Pitvd_asinh(NODE arg,Obj *rp)
1455: {
1456: mpfi_func_d(arg, mpfi_asinh, rp);
1457: }
1458:
1.6 ! kondoh 1459: static
1.5 kondoh 1460: void Pitvd_acosh(NODE arg,Obj *rp)
1461: {
1462: mpfi_func_d(arg, mpfi_acosh, rp);
1463: }
1464:
1.6 ! kondoh 1465: static
1.5 kondoh 1466: void Pitvd_atanh(NODE arg,Obj *rp)
1467: {
1468: mpfi_func_d(arg, mpfi_atanh, rp);
1469: }
1470:
1.6 ! kondoh 1471: static
1.5 kondoh 1472: void Pitvd_exp(NODE arg,Obj *rp)
1473: {
1474: mpfi_func_d(arg, mpfi_exp, rp);
1475: }
1476:
1.6 ! kondoh 1477: static
1.5 kondoh 1478: void Pitvd_log(NODE arg,Obj *rp)
1479: {
1480: mpfi_func_d(arg, mpfi_log, rp);
1481: }
1482:
1.6 ! kondoh 1483: static
1.5 kondoh 1484: void Pitvd_abs(NODE arg,Obj *rp)
1485: {
1486: mpfi_func_d(arg, mpfi_abs, rp);
1487: }
1488:
1489: /*
1490: void mp_factorial(NODE arg,Num *rp)
1491: {
1492: struct oNODE arg0;
1493: Num a,a1;
1494:
1495: a = (Num)ARG0(arg);
1496: if ( !a ) *rp = (Num)ONE;
1497: else if ( INT(a) ) Pfac(arg,rp);
1498: else {
1499: addnum(0,a,(Num)ONE,&a1);
1500: arg0.body = (pointer)a1;
1501: arg0.next = arg->next;
1502: Pmpfr_gamma(&arg0,rp);
1503: }
1504: }
1505: */
1.6 ! kondoh 1506: static
! 1507: void Pitvd_pow(NODE arg,Num *rp)
! 1508: {
! 1509: Num ii, ss;
! 1510: IntervalDouble c;
! 1511: Num rpbf;
! 1512: double inf, sup;
! 1513:
! 1514: Pitvbf_pow(arg, &rpbf);
! 1515: itvtois((Itv)rpbf, &ii, &ss);
! 1516: inf = toRealDown(ii);
! 1517: sup = toRealUp(ss);
! 1518: MKIntervalDouble(inf,sup,c);
! 1519: *rp = (Num)c;
! 1520: }
1.5 kondoh 1521:
1.6 ! kondoh 1522: static
1.5 kondoh 1523: void Pitvbf_pow(NODE arg,Num *rp)
1524: {
1525: Num a,e;
1526: BF r,re,im;
1527: C c;
1528: mpc_t mpc,a1,e1;
1529: int prec;
1530:
1531: prec = NEXT(NEXT(arg)) ? ZTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
1532: a = ARG0(arg);
1533: e = ARG1(arg);
1534: if ( !e ) {
1535: *rp = (Num)ONE;
1536: } else if ( !a ) {
1537: *rp = 0;
1538: } else if ( NID(a) == N_C || NID(e) == N_C ) {
1539: error("itv_pow() : not support args");
1540: *rp = 0;
1541: } else {
1542: Num ii, ss;
1543: Itv c;
1544: BF inf, sup;
1545: mpfi_t a_val, e_val, rv;
1546:
1547: mpfi_init2(rv,prec);
1548:
1549: itvtois((Itv)a, &ii, &ss);
1550: inf = (BF)tobf(ii, prec);
1551: sup = (BF)tobf(ss, prec);
1552: mpfi_init2(a_val,prec);
1553: mpfr_set(&(a_val->left), BDY(inf), MPFR_RNDD);
1554: mpfr_set(&(a_val->right), BDY(sup), MPFR_RNDU);
1555:
1556: itvtois((Itv)e, &ii, &ss);
1557: inf = (BF)tobf(ii, prec);
1558: sup = (BF)tobf(ss, prec);
1559: mpfi_init2(e_val,prec);
1560: mpfr_set(&(e_val->left), BDY(inf), MPFR_RNDD);
1561: mpfr_set(&(e_val->right), BDY(sup), MPFR_RNDU);
1562:
1563:
1564: mpfi_log(rv, a_val);
1565: mpfi_mul(a_val, rv, e_val);
1566: mpfi_exp(rv, a_val);
1567:
1568: MPFRTOBF(&(rv->left), inf);
1569: MPFRTOBF(&(rv->right), sup);
1570:
1571: if ( !cmpbf((Num)inf,0) ) inf = 0;
1572: if ( !cmpbf((Num)sup,0) ) sup = 0;
1573:
1574: if ( inf || sup ) {
1575: istoitv((Num)inf, (Num)sup, &c);
1576: } else {
1577: c = 0;
1578: }
1579: *rp = (Num)c;
1580: //mpfi_clear(rv);
1581: mpfi_clear(a_val);
1582: mpfi_clear(e_val);
1583: }
1584: }
1585:
1.6 ! kondoh 1586: static void Pitv_pi(NODE arg, Obj *rp)
! 1587: {
! 1588: if ( bigfloat )
! 1589: Pitvbf_pi(arg, rp);
! 1590: else
! 1591: Pitvd_pi(arg, rp);
! 1592: }
! 1593:
! 1594: static void Pitv_e(NODE arg, Obj *rp)
! 1595: {
! 1596: if ( bigfloat )
! 1597: Pitvbf_e(arg, rp);
! 1598: else
! 1599: Pitvd_e(arg, rp);
! 1600: }
! 1601:
! 1602: static void Pitv_sin(NODE arg, Obj *rp)
! 1603: {
! 1604: if ( bigfloat )
! 1605: Pitvbf_sin(arg, rp);
! 1606: else
! 1607: Pitvd_sin(arg, rp);
! 1608: }
! 1609:
! 1610: static void Pitv_cos(NODE arg, Obj *rp)
! 1611: {
! 1612: if ( bigfloat )
! 1613: Pitvbf_cos(arg, rp);
! 1614: else
! 1615: Pitvd_cos(arg, rp);
! 1616: }
! 1617:
! 1618: static void Pitv_tan(NODE arg, Obj *rp)
! 1619: {
! 1620: if ( bigfloat )
! 1621: Pitvbf_tan(arg, rp);
! 1622: else
! 1623: Pitvd_tan(arg, rp);
! 1624: }
! 1625:
! 1626: static void Pitv_asin(NODE arg, Obj *rp)
! 1627: {
! 1628: if ( bigfloat )
! 1629: Pitvbf_asin(arg, rp);
! 1630: else
! 1631: Pitvd_asin(arg, rp);
! 1632: }
! 1633:
! 1634: static void Pitv_acos(NODE arg, Obj *rp)
! 1635: {
! 1636: if ( bigfloat )
! 1637: Pitvbf_acos(arg, rp);
! 1638: else
! 1639: Pitvd_acos(arg, rp);
! 1640: }
! 1641:
! 1642: static void Pitv_atan(NODE arg, Obj *rp)
! 1643: {
! 1644: if ( bigfloat )
! 1645: Pitvbf_atan(arg, rp);
! 1646: else
! 1647: Pitvd_atan(arg, rp);
! 1648: }
! 1649:
! 1650: static void Pitv_sinh(NODE arg, Obj *rp)
! 1651: {
! 1652: if ( bigfloat )
! 1653: Pitvbf_sinh(arg, rp);
! 1654: else
! 1655: Pitvd_sinh(arg, rp);
! 1656: }
! 1657:
! 1658: static void Pitv_cosh(NODE arg, Obj *rp)
! 1659: {
! 1660: if ( bigfloat )
! 1661: Pitvbf_cosh(arg, rp);
! 1662: else
! 1663: Pitvd_cosh(arg, rp);
! 1664: }
! 1665:
! 1666: static void Pitv_tanh(NODE arg, Obj *rp)
! 1667: {
! 1668: if ( bigfloat )
! 1669: Pitvbf_tanh(arg, rp);
! 1670: else
! 1671: Pitvd_tanh(arg, rp);
! 1672: }
! 1673:
! 1674: static void Pitv_asinh(NODE arg, Obj *rp)
! 1675: {
! 1676: if ( bigfloat )
! 1677: Pitvbf_asinh(arg, rp);
! 1678: else
! 1679: Pitvd_asinh(arg, rp);
! 1680: }
! 1681:
! 1682: static void Pitv_acosh(NODE arg, Obj *rp)
! 1683: {
! 1684: if ( bigfloat )
! 1685: Pitvbf_acosh(arg, rp);
! 1686: else
! 1687: Pitvd_acosh(arg, rp);
! 1688: }
! 1689:
! 1690: static void Pitv_atanh(NODE arg, Obj *rp)
! 1691: {
! 1692: if ( bigfloat )
! 1693: Pitvbf_atanh(arg, rp);
! 1694: else
! 1695: Pitvd_atanh(arg, rp);
! 1696: }
! 1697:
! 1698: static void Pitv_exp(NODE arg, Obj *rp)
! 1699: {
! 1700: if ( bigfloat )
! 1701: Pitvbf_exp(arg, rp);
! 1702: else
! 1703: Pitvd_exp(arg, rp);
! 1704: }
! 1705:
! 1706: static void Pitv_log(NODE arg, Obj *rp)
! 1707: {
! 1708: if ( bigfloat )
! 1709: Pitvbf_log(arg, rp);
! 1710: else
! 1711: Pitvd_log(arg, rp);
! 1712: }
! 1713:
! 1714: static void Pitv_abs(NODE arg, Obj *rp)
! 1715: {
! 1716: if ( bigfloat )
! 1717: Pitvbf_abs(arg, rp);
! 1718: else
! 1719: Pitvd_abs(arg, rp);
! 1720: }
! 1721:
! 1722: static void Pitv_pow(NODE arg, Num *rp)
! 1723: {
! 1724: if ( bigfloat )
! 1725: Pitvbf_pow(arg, rp);
! 1726: else
! 1727: Pitvd_pow(arg, rp);
! 1728: }
! 1729:
1.5 kondoh 1730: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to itvnum.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / builtin