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