Annotation of OpenXM_contrib2/asir2000/builtin/bfaux.c, Revision 1.10
1.10 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.9 2015/08/17 05:18:36 noro Exp $ */
1.1 noro 2: #include "ca.h"
3: #include "parse.h"
4:
1.4 noro 5: void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround();
1.8 noro 6: void Pmpfr_ai();
1.9 noro 7: void Pmpfr_eint(), Pmpfr_erf(), Pmpfr_erfc(), Pmpfr_li2();
1.8 noro 8: void Pmpfr_zeta();
9: void Pmpfr_j0(), Pmpfr_j1();
10: void Pmpfr_y0(), Pmpfr_y1();
11: void Pmpfr_gamma(), Pmpfr_lngamma(), Pmpfr_digamma();
12: void Pmpfr_floor(), Pmpfr_round(), Pmpfr_ceil();
1.1 noro 13:
14: struct ftab bf_tab[] = {
15: {"eval",Peval,-2},
16: {"setprec",Psetprec,-1},
1.3 noro 17: {"setbprec",Psetbprec,-1},
1.4 noro 18: {"setround",Psetround,-1},
1.1 noro 19: {"todouble",Ptodouble,1},
1.8 noro 20: {"mpfr_ai",Pmpfr_ai,-2},
21: {"mpfr_zeta",Pmpfr_zeta,-2},
22: {"mpfr_j0",Pmpfr_j0,-2},
23: {"mpfr_j1",Pmpfr_j1,-2},
24: {"mpfr_y0",Pmpfr_y0,-2},
25: {"mpfr_y1",Pmpfr_y1,-2},
26: {"mpfr_eint",Pmpfr_eint,-2},
27: {"mpfr_erf",Pmpfr_erf,-2},
1.9 noro 28: {"mpfr_erfc",Pmpfr_erfc,-2},
1.8 noro 29: {"mpfr_li2",Pmpfr_li2,-2},
1.5 noro 30: {"mpfr_gamma",Pmpfr_gamma,-2},
1.8 noro 31: {"mpfr_lngamma",Pmpfr_gamma,-2},
32: {"mpfr_digamma",Pmpfr_gamma,-2},
33: {"mpfr_floor",Pmpfr_floor,-2},
34: {"mpfr_ceil",Pmpfr_ceil,-2},
35: {"mpfr_round",Pmpfr_round,-2},
1.1 noro 36: {0,0,0},
37: };
38:
1.4 noro 39: int mpfr_roundmode = MPFR_RNDN;
40:
1.1 noro 41: void Ptodouble(NODE arg,Num *rp)
42: {
43: double r,i;
44: Real real,imag;
45: Num num;
46:
47: asir_assert(ARG0(arg),O_N,"todouble");
48: num = (Num)ARG0(arg);
49: if ( !num ) {
50: *rp = 0;
51: return;
52: }
53: switch ( NID(num) ) {
54: case N_R: case N_Q: case N_B:
55: r = ToReal(num);
56: MKReal(r,real);
57: *rp = (Num)real;
58: break;
59: case N_C:
60: r = ToReal(((C)num)->r);
61: i = ToReal(((C)num)->i);
62: MKReal(r,real);
63: MKReal(i,imag);
64: reimtocplx((Num)real,(Num)imag,rp);
65: break;
66: default:
67: *rp = num;
68: break;
69: }
70: }
71:
1.4 noro 72: void Peval(NODE arg,Obj *rp)
1.1 noro 73: {
74: int prec;
75:
76: asir_assert(ARG0(arg),O_R,"eval");
77: if ( argc(arg) == 2 ) {
1.3 noro 78: prec = QTOS((Q)ARG1(arg))*3.32193;
1.1 noro 79: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
80: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
81: } else
82: prec = 0;
1.2 noro 83: evalr(CO,(Obj)ARG0(arg),prec,rp);
1.1 noro 84: }
85:
1.3 noro 86: /* set/get decimal precision */
1.1 noro 87:
88: void Psetprec(NODE arg,Obj *rp)
89: {
90: int p;
91: Q q;
1.3 noro 92: int prec,dprec;
93:
94: prec = mpfr_get_default_prec();
95: /* decimal precision */
96: dprec = prec*0.30103;
97: STOQ(dprec,q); *rp = (Obj)q;
98: if ( arg ) {
99: asir_assert(ARG0(arg),O_N,"setprec");
100: prec = QTOS((Q)ARG0(arg))*3.32193;
101: if ( p > 0 )
102: prec = p;
103: }
104: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
105: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
106: mpfr_set_default_prec(prec);
107: }
1.1 noro 108:
1.3 noro 109: /* set/get bit precision */
1.1 noro 110:
1.3 noro 111: void Psetbprec(NODE arg,Obj *rp)
112: {
113: int p;
114: Q q;
115: int prec;
116:
117: prec = mpfr_get_default_prec();
118: STOQ(prec,q); *rp = (Obj)q;
1.1 noro 119: if ( arg ) {
1.3 noro 120: asir_assert(ARG0(arg),O_N,"setbprec");
121: prec = QTOS((Q)ARG0(arg));
1.1 noro 122: if ( p > 0 )
123: prec = p;
124: }
125: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
126: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
127: mpfr_set_default_prec(prec);
128: }
129:
1.4 noro 130: void Psetround(NODE arg,Q *rp)
131: {
132: int round;
133:
134: STOQ(mpfr_roundmode,*rp);
135: if ( arg ) {
136: asir_assert(ARG0(arg),O_N,"setround");
137: round = QTOS((Q)ARG0(arg));
138: switch ( round ) {
139: case 0:
140: mpfr_roundmode = MPFR_RNDN;
141: break;
142: case 1:
143: mpfr_roundmode = MPFR_RNDZ;
144: break;
145: case 2:
146: mpfr_roundmode = MPFR_RNDU;
147: break;
148: case 3:
149: mpfr_roundmode = MPFR_RNDD;
150: break;
151: case 4:
152: mpfr_roundmode = MPFR_RNDA;
153: break;
154: case 5:
155: mpfr_roundmode = MPFR_RNDF;
156: break;
157: case 6:
158: mpfr_roundmode = MPFR_RNDNA;
159: break;
160: default:
161: error("setround : invalid rounding mode");
162: break;
163: }
164: }
165: }
166:
1.1 noro 167: Num tobf(Num a,int prec);
168:
169: void mp_pi(NODE arg,BF *rp)
170: {
1.10 ! noro 171: int prec;
1.1 noro 172: BF r;
173:
174: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
175: NEWBF(r);
176: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 177: mpfr_const_pi(r->body,mpfr_roundmode);
1.10 ! noro 178: if ( !cmpbf((Num)r,0) ) r = 0;
! 179: *rp = r;
1.1 noro 180: }
181:
182: void mp_e(NODE arg,BF *rp)
183: {
1.10 ! noro 184: int prec;
1.1 noro 185: mpfr_t one;
186: BF r;
187:
188: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
189: NEWBF(r);
190: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
191: mpfr_init(one);
1.4 noro 192: mpfr_set_ui(one,1,mpfr_roundmode);
193: mpfr_exp(r->body,one,mpfr_roundmode);
1.10 ! noro 194: if ( !cmpbf((Num)r,0) ) r = 0;
! 195: *rp = r;
1.1 noro 196: }
197:
1.10 ! noro 198: void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)(),Num *rp)
1.1 noro 199: {
200: Num a;
1.10 ! noro 201: int prec;
! 202: BF r,re,im;
! 203: C c;
! 204: mpc_t mpc,a1;
1.1 noro 205:
1.10 ! noro 206: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
1.1 noro 207: a = tobf(ARG0(arg),prec);
1.10 ! noro 208: if ( a && NID(a)==N_C ) {
! 209: mpc_init2(mpc,prec); mpc_init2(a1,prec);
! 210: re = (BF)((C)a)->r; im = (BF)((C)a)->i;
! 211: mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
! 212: (*mpc_f)(mpc,a1,mpfr_roundmode);
! 213: MPFRTOBF(mpc_realref(mpc),re);
! 214: MPFRTOBF(mpc_imagref(mpc),im);
! 215: if ( !cmpbf((Num)re,0) ) re = 0;
! 216: if ( !cmpbf((Num)im,0) ) im = 0;
! 217: if ( !im )
! 218: *rp = (Num)re;
! 219: else {
! 220: NEWC(c); c->r = (Num)re; c->i = (Num)im;
! 221: *rp = (Num)c;
! 222: }
! 223: } else {
! 224: NEWBF(r);
! 225: mpfr_init2(r->body,prec);
! 226: (*mpfr_f)(r->body,((BF)a)->body,mpfr_roundmode);
! 227: if ( !cmpbf((Num)r,0) ) r = 0;
! 228: *rp = (Num)r;
! 229: }
1.1 noro 230: }
231:
1.10 ! noro 232: void mp_sin(NODE arg,Num *rp)
1.1 noro 233: {
1.10 ! noro 234: mpfr_or_mpc(arg,mpfr_sin,mpc_sin,rp);
! 235: }
1.1 noro 236:
1.10 ! noro 237: void mp_cos(NODE arg,Num *rp)
! 238: {
! 239: mpfr_or_mpc(arg,mpfr_cos,mpc_cos,rp);
1.1 noro 240: }
241:
1.10 ! noro 242: void mp_tan(NODE arg,Num *rp)
1.1 noro 243: {
1.10 ! noro 244: mpfr_or_mpc(arg,mpfr_tan,mpc_tan,rp);
! 245: }
1.1 noro 246:
1.10 ! noro 247: void mp_asin(NODE arg,Num *rp)
! 248: {
! 249: mpfr_or_mpc(arg,mpfr_asin,mpc_asin,rp);
1.1 noro 250: }
251:
1.10 ! noro 252: void mp_acos(NODE arg,Num *rp)
1.1 noro 253: {
1.10 ! noro 254: mpfr_or_mpc(arg,mpfr_acos,mpc_acos,rp);
! 255: }
1.1 noro 256:
1.10 ! noro 257: void mp_atan(NODE arg,Num *rp)
! 258: {
! 259: mpfr_or_mpc(arg,mpfr_atan,mpc_atan,rp);
1.1 noro 260: }
261:
1.10 ! noro 262: void mp_sinh(NODE arg,Num *rp)
1.1 noro 263: {
1.10 ! noro 264: mpfr_or_mpc(arg,mpfr_sinh,mpc_sinh,rp);
1.1 noro 265: }
266:
1.10 ! noro 267: void mp_cosh(NODE arg,Num *rp)
1.1 noro 268: {
1.10 ! noro 269: mpfr_or_mpc(arg,mpfr_cosh,mpc_cosh,rp);
1.1 noro 270: }
271:
1.10 ! noro 272: void mp_tanh(NODE arg,Num *rp)
1.1 noro 273: {
1.10 ! noro 274: mpfr_or_mpc(arg,mpfr_tanh,mpc_tanh,rp);
1.1 noro 275: }
276:
1.10 ! noro 277: void mp_asinh(NODE arg,Num *rp)
1.1 noro 278: {
1.10 ! noro 279: mpfr_or_mpc(arg,mpfr_asinh,mpc_asinh,rp);
1.1 noro 280: }
281:
1.10 ! noro 282: void mp_acosh(NODE arg,Num *rp)
1.1 noro 283: {
1.10 ! noro 284: mpfr_or_mpc(arg,mpfr_acosh,mpc_acosh,rp);
1.1 noro 285: }
286:
1.10 ! noro 287: void mp_atanh(NODE arg,Num *rp)
1.1 noro 288: {
1.10 ! noro 289: mpfr_or_mpc(arg,mpfr_atanh,mpc_atanh,rp);
1.1 noro 290: }
291:
1.10 ! noro 292: void mp_exp(NODE arg,Num *rp)
1.1 noro 293: {
1.10 ! noro 294: mpfr_or_mpc(arg,mpfr_exp,mpc_exp,rp);
1.1 noro 295: }
296:
1.10 ! noro 297: void mp_log(NODE arg,Num *rp)
1.1 noro 298: {
1.10 ! noro 299: mpfr_or_mpc(arg,mpfr_log,mpc_log,rp);
1.1 noro 300: }
301:
1.10 ! noro 302: void mp_pow(NODE arg,Num *rp)
1.1 noro 303: {
304: Num a,e;
1.10 ! noro 305: int prec;
! 306: BF r,re,im;
! 307: C c;
! 308: mpc_t mpc,a1,e1;
1.1 noro 309:
1.10 ! noro 310: prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
1.1 noro 311: a = tobf(ARG0(arg),prec);
312: e = tobf(ARG1(arg),prec);
1.10 ! noro 313: if ( NID(a) == N_C || NID(e) == N_C || MPFR_SIGN(((BF)a)->body) < 0 ) {
! 314: mpc_init2(mpc,prec); mpc_init2(a1,prec); mpc_init2(e1,prec);
! 315: if ( NID(a) == N_C ) {
! 316: re = (BF)((C)a)->r; im = (BF)((C)a)->i;
! 317: mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
! 318: } else {
! 319: re = (BF)a;
! 320: mpc_set_fr(a1,re->body,mpfr_roundmode);
! 321: }
! 322: if ( NID(e) == N_C ) {
! 323: re = (BF)((C)e)->r; im = (BF)((C)e)->i;
! 324: mpc_set_fr_fr(e1,re->body,im->body,mpfr_roundmode);
! 325: } else {
! 326: re = (BF)e;
! 327: mpc_set_fr(e1,re->body,mpfr_roundmode);
! 328: }
! 329: mpc_pow(mpc,a1,e1,mpfr_roundmode);
! 330: MPFRTOBF(mpc_realref(mpc),re);
! 331: MPFRTOBF(mpc_imagref(mpc),im);
! 332: if ( !cmpbf((Num)re,0) ) re = 0;
! 333: if ( !cmpbf((Num)im,0) ) im = 0;
! 334: if ( !im )
! 335: *rp = (Num)re;
! 336: else {
! 337: NEWC(c); c->r = (Num)re; c->i = (Num)im;
! 338: *rp = (Num)c;
! 339: }
! 340: } else {
! 341: NEWBF(r);
! 342: mpfr_init2(r->body,prec);
! 343: mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
! 344: *rp = (Num)r;
! 345: }
1.1 noro 346: }
1.5 noro 347:
1.8 noro 348: #define SETPREC \
349: (prec)=NEXT(arg)?QTOS((Q)ARG1(arg)):0;\
350: (prec)*=3.32193;\
351: (a)=tobf(ARG0(arg),prec);\
352: NEWBF(r);\
353: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
354:
355:
1.5 noro 356: void Pmpfr_gamma(NODE arg,BF *rp)
357: {
358: Num a;
359: int prec;
360: BF r;
361:
1.8 noro 362: SETPREC
1.5 noro 363: mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
364: *rp = r;
365: }
1.6 takayama 366:
1.8 noro 367: void Pmpfr_lngamma(NODE arg,BF *rp)
368: {
369: Num a;
370: int prec;
371: BF r;
372:
373: SETPREC
374: mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode);
375: *rp = r;
376: }
377:
378: void Pmpfr_digamma(NODE arg,BF *rp)
379: {
380: Num a;
381: int prec;
382: BF r;
383:
384: SETPREC
385: mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode);
386: *rp = r;
387: }
388:
389: void Pmpfr_zeta(NODE arg,BF *rp)
390: {
391: Num a;
392: int prec;
393: BF r;
394:
395: SETPREC
396: mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode);
397: *rp = r;
398: }
399:
400: void Pmpfr_eint(NODE arg,BF *rp)
401: {
402: Num a;
403: int prec;
404: BF r;
405:
406: SETPREC
407: mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode);
408: *rp = r;
409: }
410:
411: void Pmpfr_erf(NODE arg,BF *rp)
412: {
413: Num a;
414: int prec;
415: BF r;
416:
417: SETPREC
418: mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode);
419: *rp = r;
420: }
421:
1.9 noro 422: void Pmpfr_erfc(NODE arg,BF *rp)
423: {
424: Num a;
425: int prec;
426: BF r;
427:
428: SETPREC
429: mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode);
430: *rp = r;
431: }
432:
1.8 noro 433: void Pmpfr_j0(NODE arg,BF *rp)
434: {
435: Num a;
436: int prec;
437: BF r;
438:
439: SETPREC
440: mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode);
441: *rp = r;
442: }
443:
444: void Pmpfr_j1(NODE arg,BF *rp)
445: {
446: Num a;
447: int prec;
448: BF r;
449:
450: SETPREC
451: mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode);
452: *rp = r;
453: }
454:
455: void Pmpfr_y0(NODE arg,BF *rp)
456: {
457: Num a;
458: int prec;
459: BF r;
460:
461: SETPREC
462: mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode);
463: *rp = r;
464: }
465:
466: void Pmpfr_y1(NODE arg,BF *rp)
467: {
468: Num a;
469: int prec;
470: BF r;
471:
472: SETPREC
473: mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode);
474: *rp = r;
475: }
476:
477: void Pmpfr_li2(NODE arg,BF *rp)
478: {
479: Num a;
480: int prec;
481: BF r;
482:
483: SETPREC
484: mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode);
485: *rp = r;
486: }
487:
488: void Pmpfr_ai(NODE arg,BF *rp)
489: {
490: Num a;
491: int prec;
492: BF r;
493:
494: SETPREC
495: mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode);
496: *rp = r;
497: }
498:
1.7 takayama 499: void Pmpfr_floor(NODE arg,Q *rp)
1.6 takayama 500: {
501: Num a;
502: int prec;
503: BF r;
1.7 takayama 504: mpz_t t;
505: GZ rz;
1.6 takayama 506:
1.8 noro 507: SETPREC
1.6 takayama 508: mpfr_floor(r->body,((BF)a)->body);
1.7 takayama 509: mpz_init(t);
510: mpfr_get_z(t,r->body,mpfr_roundmode);
511: MPZTOGZ(t,rz);
512: *rp = gztoz(rz);
513: }
514:
1.8 noro 515: void Pmpfr_ceil(NODE arg,Q *rp)
516: {
517: Num a;
518: int prec;
519: BF r;
520: mpz_t t;
521: GZ rz;
522:
523: SETPREC
524: mpfr_ceil(r->body,((BF)a)->body);
525: mpz_init(t);
526: mpfr_get_z(t,r->body,mpfr_roundmode);
527: MPZTOGZ(t,rz);
528: *rp = gztoz(rz);
529: }
530:
1.7 takayama 531: void Pmpfr_round(NODE arg,Q *rp)
532: {
533: Num a;
534: int prec;
535: BF r;
536: mpz_t t;
537: GZ rz;
538:
1.8 noro 539: SETPREC
1.7 takayama 540: mpfr_round(r->body,((BF)a)->body);
541: mpz_init(t);
542: mpfr_get_z(t,r->body,mpfr_roundmode);
543: MPZTOGZ(t,rz);
544: *rp = gztoz(rz);
1.6 takayama 545: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bfaux.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin