![[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|
Return to bfaux.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.9 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.8 2015/08/16 03:12:09 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: {
171: int prec;
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.1 noro 178: *rp = r;
179: }
180:
181: void mp_e(NODE arg,BF *rp)
182: {
183: int prec;
184: mpfr_t one;
185: BF r;
186:
187: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
188: NEWBF(r);
189: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
190: mpfr_init(one);
1.4 noro 191: mpfr_set_ui(one,1,mpfr_roundmode);
192: mpfr_exp(r->body,one,mpfr_roundmode);
1.1 noro 193: *rp = r;
194: }
195:
196: void mp_sin(NODE arg,BF *rp)
197: {
198: Num a;
199: int prec;
200: BF r;
201:
202: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
203: a = tobf(ARG0(arg),prec);
204: NEWBF(r);
205: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 206: mpfr_sin(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 207: *rp = r;
208: }
209:
210: void mp_cos(NODE arg,BF *rp)
211: {
212: Num a;
213: int prec;
214: BF r;
215:
216: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
217: a = tobf(ARG0(arg),prec);
218: NEWBF(r);
219: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 220: mpfr_cos(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 221: *rp = r;
222: }
223:
224: void mp_tan(NODE arg,BF *rp)
225: {
226: Num a;
227: int prec;
228: BF r;
229:
230: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
231: a = tobf(ARG0(arg),prec);
232: NEWBF(r);
233: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 234: mpfr_tan(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 235: *rp = r;
236: }
237:
238: void mp_asin(NODE arg,BF *rp)
239: {
240: Num a;
241: int prec;
242: BF r;
243:
244: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
245: a = tobf(ARG0(arg),prec);
246: NEWBF(r);
247: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 248: mpfr_asin(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 249: *rp = r;
250: }
251: void mp_acos(NODE arg,BF *rp)
252: {
253: Num a;
254: int prec;
255: BF r;
256:
257: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
258: a = tobf(ARG0(arg),prec);
259: NEWBF(r);
260: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 261: mpfr_acos(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 262: *rp = r;
263: }
264: void mp_atan(NODE arg,BF *rp)
265: {
266: Num a;
267: int prec;
268: BF r;
269:
270: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
271: a = tobf(ARG0(arg),prec);
272: NEWBF(r);
273: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 274: mpfr_atan(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 275: *rp = r;
276: }
277:
278: void mp_sinh(NODE arg,BF *rp)
279: {
280: Num a;
281: int prec;
282: BF r;
283:
284: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
285: a = tobf(ARG0(arg),prec);
286: NEWBF(r);
287: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 288: mpfr_sinh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 289: *rp = r;
290: }
291:
292: void mp_cosh(NODE arg,BF *rp)
293: {
294: Num a;
295: int prec;
296: BF r;
297:
298: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
299: a = tobf(ARG0(arg),prec);
300: NEWBF(r);
301: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 302: mpfr_cosh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 303: *rp = r;
304: }
305:
306: void mp_tanh(NODE arg,BF *rp)
307: {
308: Num a;
309: int prec;
310: BF r;
311:
312: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
313: a = tobf(ARG0(arg),prec);
314: NEWBF(r);
315: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 316: mpfr_tanh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 317: *rp = r;
318: }
319:
320: void mp_asinh(NODE arg,BF *rp)
321: {
322: Num a;
323: int prec;
324: BF r;
325:
326: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
327: a = tobf(ARG0(arg),prec);
328: NEWBF(r);
329: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 330: mpfr_asinh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 331: *rp = r;
332: }
333: void mp_acosh(NODE arg,BF *rp)
334: {
335: Num a;
336: int prec;
337: BF r;
338:
339: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
340: a = tobf(ARG0(arg),prec);
341: NEWBF(r);
342: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 343: mpfr_acosh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 344: *rp = r;
345: }
346: void mp_atanh(NODE arg,BF *rp)
347: {
348: Num a;
349: int prec;
350: BF r;
351:
352: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
353: a = tobf(ARG0(arg),prec);
354: NEWBF(r);
355: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 356: mpfr_atanh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 357: *rp = r;
358: }
359:
360: void mp_exp(NODE arg,BF *rp)
361: {
362: Num a;
363: int prec;
364: BF r;
365:
366: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
367: a = tobf(ARG0(arg),prec);
368: NEWBF(r);
369: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 370: mpfr_exp(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 371: *rp = r;
372: }
373:
374: void mp_log(NODE arg,BF *rp)
375: {
376: Num a;
377: int prec;
378: BF r;
379:
380: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
381: a = tobf(ARG0(arg),prec);
382: NEWBF(r);
383: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 384: mpfr_log(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 385: *rp = r;
386: }
387:
388: void mp_pow(NODE arg,BF *rp)
389: {
390: Num a,e;
391: int prec;
392: BF r;
393:
394: prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : 0;
395: a = tobf(ARG0(arg),prec);
396: e = tobf(ARG1(arg),prec);
397: NEWBF(r);
398: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 399: mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
1.1 noro 400: *rp = r;
401: }
1.5 noro 402:
1.8 noro 403: #define SETPREC \
404: (prec)=NEXT(arg)?QTOS((Q)ARG1(arg)):0;\
405: (prec)*=3.32193;\
406: (a)=tobf(ARG0(arg),prec);\
407: NEWBF(r);\
408: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
409:
410:
1.5 noro 411: void Pmpfr_gamma(NODE arg,BF *rp)
412: {
413: Num a;
414: int prec;
415: BF r;
416:
1.8 noro 417: SETPREC
1.5 noro 418: mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
419: *rp = r;
420: }
1.6 takayama 421:
1.8 noro 422: void Pmpfr_lngamma(NODE arg,BF *rp)
423: {
424: Num a;
425: int prec;
426: BF r;
427:
428: SETPREC
429: mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode);
430: *rp = r;
431: }
432:
433: void Pmpfr_digamma(NODE arg,BF *rp)
434: {
435: Num a;
436: int prec;
437: BF r;
438:
439: SETPREC
440: mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode);
441: *rp = r;
442: }
443:
444: void Pmpfr_zeta(NODE arg,BF *rp)
445: {
446: Num a;
447: int prec;
448: BF r;
449:
450: SETPREC
451: mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode);
452: *rp = r;
453: }
454:
455: void Pmpfr_eint(NODE arg,BF *rp)
456: {
457: Num a;
458: int prec;
459: BF r;
460:
461: SETPREC
462: mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode);
463: *rp = r;
464: }
465:
466: void Pmpfr_erf(NODE arg,BF *rp)
467: {
468: Num a;
469: int prec;
470: BF r;
471:
472: SETPREC
473: mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode);
474: *rp = r;
475: }
476:
1.9 ! noro 477: void Pmpfr_erfc(NODE arg,BF *rp)
! 478: {
! 479: Num a;
! 480: int prec;
! 481: BF r;
! 482:
! 483: SETPREC
! 484: mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode);
! 485: *rp = r;
! 486: }
! 487:
1.8 noro 488: void Pmpfr_j0(NODE arg,BF *rp)
489: {
490: Num a;
491: int prec;
492: BF r;
493:
494: SETPREC
495: mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode);
496: *rp = r;
497: }
498:
499: void Pmpfr_j1(NODE arg,BF *rp)
500: {
501: Num a;
502: int prec;
503: BF r;
504:
505: SETPREC
506: mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode);
507: *rp = r;
508: }
509:
510: void Pmpfr_y0(NODE arg,BF *rp)
511: {
512: Num a;
513: int prec;
514: BF r;
515:
516: SETPREC
517: mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode);
518: *rp = r;
519: }
520:
521: void Pmpfr_y1(NODE arg,BF *rp)
522: {
523: Num a;
524: int prec;
525: BF r;
526:
527: SETPREC
528: mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode);
529: *rp = r;
530: }
531:
532: void Pmpfr_li2(NODE arg,BF *rp)
533: {
534: Num a;
535: int prec;
536: BF r;
537:
538: SETPREC
539: mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode);
540: *rp = r;
541: }
542:
543: void Pmpfr_ai(NODE arg,BF *rp)
544: {
545: Num a;
546: int prec;
547: BF r;
548:
549: SETPREC
550: mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode);
551: *rp = r;
552: }
553:
1.7 takayama 554: void Pmpfr_floor(NODE arg,Q *rp)
1.6 takayama 555: {
556: Num a;
557: int prec;
558: BF r;
1.7 takayama 559: mpz_t t;
560: GZ rz;
1.6 takayama 561:
1.8 noro 562: SETPREC
1.6 takayama 563: mpfr_floor(r->body,((BF)a)->body);
1.7 takayama 564: mpz_init(t);
565: mpfr_get_z(t,r->body,mpfr_roundmode);
566: MPZTOGZ(t,rz);
567: *rp = gztoz(rz);
568: }
569:
1.8 noro 570: void Pmpfr_ceil(NODE arg,Q *rp)
571: {
572: Num a;
573: int prec;
574: BF r;
575: mpz_t t;
576: GZ rz;
577:
578: SETPREC
579: mpfr_ceil(r->body,((BF)a)->body);
580: mpz_init(t);
581: mpfr_get_z(t,r->body,mpfr_roundmode);
582: MPZTOGZ(t,rz);
583: *rp = gztoz(rz);
584: }
585:
1.7 takayama 586: void Pmpfr_round(NODE arg,Q *rp)
587: {
588: Num a;
589: int prec;
590: BF r;
591: mpz_t t;
592: GZ rz;
593:
1.8 noro 594: SETPREC
1.7 takayama 595: mpfr_round(r->body,((BF)a)->body);
596: mpz_init(t);
597: mpfr_get_z(t,r->body,mpfr_roundmode);
598: MPZTOGZ(t,rz);
599: *rp = gztoz(rz);
1.6 takayama 600: }