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