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