Annotation of OpenXM_contrib2/asir2000/builtin/bfaux.c, Revision 1.7
1.7 ! takayama 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.6 2015/08/07 06:15:00 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.7 ! takayama 6: void Pmpfr_gamma(), Pmpfr_floor(), Pmpfr_round();
1.1 noro 7:
8: struct ftab bf_tab[] = {
9: {"eval",Peval,-2},
10: {"setprec",Psetprec,-1},
1.3 noro 11: {"setbprec",Psetbprec,-1},
1.4 noro 12: {"setround",Psetround,-1},
1.1 noro 13: {"todouble",Ptodouble,1},
1.5 noro 14: {"mpfr_gamma",Pmpfr_gamma,-2},
1.6 takayama 15: {"mpfr_floor",Pmpfr_floor,-1},
1.7 ! takayama 16: {"mpfr_round",Pmpfr_round,-1},
1.1 noro 17: {0,0,0},
18: };
19:
1.4 noro 20: int mpfr_roundmode = MPFR_RNDN;
21:
1.1 noro 22: void Ptodouble(NODE arg,Num *rp)
23: {
24: double r,i;
25: Real real,imag;
26: Num num;
27:
28: asir_assert(ARG0(arg),O_N,"todouble");
29: num = (Num)ARG0(arg);
30: if ( !num ) {
31: *rp = 0;
32: return;
33: }
34: switch ( NID(num) ) {
35: case N_R: case N_Q: case N_B:
36: r = ToReal(num);
37: MKReal(r,real);
38: *rp = (Num)real;
39: break;
40: case N_C:
41: r = ToReal(((C)num)->r);
42: i = ToReal(((C)num)->i);
43: MKReal(r,real);
44: MKReal(i,imag);
45: reimtocplx((Num)real,(Num)imag,rp);
46: break;
47: default:
48: *rp = num;
49: break;
50: }
51: }
52:
1.4 noro 53: void Peval(NODE arg,Obj *rp)
1.1 noro 54: {
55: int prec;
56:
57: asir_assert(ARG0(arg),O_R,"eval");
58: if ( argc(arg) == 2 ) {
1.3 noro 59: prec = QTOS((Q)ARG1(arg))*3.32193;
1.1 noro 60: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
61: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
62: } else
63: prec = 0;
1.2 noro 64: evalr(CO,(Obj)ARG0(arg),prec,rp);
1.1 noro 65: }
66:
1.3 noro 67: /* set/get decimal precision */
1.1 noro 68:
69: void Psetprec(NODE arg,Obj *rp)
70: {
71: int p;
72: Q q;
1.3 noro 73: int prec,dprec;
74:
75: prec = mpfr_get_default_prec();
76: /* decimal precision */
77: dprec = prec*0.30103;
78: STOQ(dprec,q); *rp = (Obj)q;
79: if ( arg ) {
80: asir_assert(ARG0(arg),O_N,"setprec");
81: prec = QTOS((Q)ARG0(arg))*3.32193;
82: if ( p > 0 )
83: prec = p;
84: }
85: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
86: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
87: mpfr_set_default_prec(prec);
88: }
1.1 noro 89:
1.3 noro 90: /* set/get bit precision */
1.1 noro 91:
1.3 noro 92: void Psetbprec(NODE arg,Obj *rp)
93: {
94: int p;
95: Q q;
96: int prec;
97:
98: prec = mpfr_get_default_prec();
99: STOQ(prec,q); *rp = (Obj)q;
1.1 noro 100: if ( arg ) {
1.3 noro 101: asir_assert(ARG0(arg),O_N,"setbprec");
102: prec = QTOS((Q)ARG0(arg));
1.1 noro 103: if ( p > 0 )
104: prec = p;
105: }
106: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
107: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
108: mpfr_set_default_prec(prec);
109: }
110:
1.4 noro 111: void Psetround(NODE arg,Q *rp)
112: {
113: int round;
114:
115: STOQ(mpfr_roundmode,*rp);
116: if ( arg ) {
117: asir_assert(ARG0(arg),O_N,"setround");
118: round = QTOS((Q)ARG0(arg));
119: switch ( round ) {
120: case 0:
121: mpfr_roundmode = MPFR_RNDN;
122: break;
123: case 1:
124: mpfr_roundmode = MPFR_RNDZ;
125: break;
126: case 2:
127: mpfr_roundmode = MPFR_RNDU;
128: break;
129: case 3:
130: mpfr_roundmode = MPFR_RNDD;
131: break;
132: case 4:
133: mpfr_roundmode = MPFR_RNDA;
134: break;
135: case 5:
136: mpfr_roundmode = MPFR_RNDF;
137: break;
138: case 6:
139: mpfr_roundmode = MPFR_RNDNA;
140: break;
141: default:
142: error("setround : invalid rounding mode");
143: break;
144: }
145: }
146: }
147:
1.1 noro 148: Num tobf(Num a,int prec);
149:
150: void mp_pi(NODE arg,BF *rp)
151: {
152: int prec;
153: BF r;
154:
155: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
156: NEWBF(r);
157: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 158: mpfr_const_pi(r->body,mpfr_roundmode);
1.1 noro 159: *rp = r;
160: }
161:
162: void mp_e(NODE arg,BF *rp)
163: {
164: int prec;
165: mpfr_t one;
166: BF r;
167:
168: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
169: NEWBF(r);
170: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
171: mpfr_init(one);
1.4 noro 172: mpfr_set_ui(one,1,mpfr_roundmode);
173: mpfr_exp(r->body,one,mpfr_roundmode);
1.1 noro 174: *rp = r;
175: }
176:
177: void mp_sin(NODE arg,BF *rp)
178: {
179: Num a;
180: int prec;
181: BF r;
182:
183: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
184: a = tobf(ARG0(arg),prec);
185: NEWBF(r);
186: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 187: mpfr_sin(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 188: *rp = r;
189: }
190:
191: void mp_cos(NODE arg,BF *rp)
192: {
193: Num a;
194: int prec;
195: BF r;
196:
197: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
198: a = tobf(ARG0(arg),prec);
199: NEWBF(r);
200: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 201: mpfr_cos(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 202: *rp = r;
203: }
204:
205: void mp_tan(NODE arg,BF *rp)
206: {
207: Num a;
208: int prec;
209: BF r;
210:
211: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
212: a = tobf(ARG0(arg),prec);
213: NEWBF(r);
214: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 215: mpfr_tan(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 216: *rp = r;
217: }
218:
219: void mp_asin(NODE arg,BF *rp)
220: {
221: Num a;
222: int prec;
223: BF r;
224:
225: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
226: a = tobf(ARG0(arg),prec);
227: NEWBF(r);
228: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 229: mpfr_asin(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 230: *rp = r;
231: }
232: void mp_acos(NODE arg,BF *rp)
233: {
234: Num a;
235: int prec;
236: BF r;
237:
238: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
239: a = tobf(ARG0(arg),prec);
240: NEWBF(r);
241: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 242: mpfr_acos(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 243: *rp = r;
244: }
245: void mp_atan(NODE arg,BF *rp)
246: {
247: Num a;
248: int prec;
249: BF r;
250:
251: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
252: a = tobf(ARG0(arg),prec);
253: NEWBF(r);
254: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 255: mpfr_atan(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 256: *rp = r;
257: }
258:
259: void mp_sinh(NODE arg,BF *rp)
260: {
261: Num a;
262: int prec;
263: BF r;
264:
265: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
266: a = tobf(ARG0(arg),prec);
267: NEWBF(r);
268: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 269: mpfr_sinh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 270: *rp = r;
271: }
272:
273: void mp_cosh(NODE arg,BF *rp)
274: {
275: Num a;
276: int prec;
277: BF r;
278:
279: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
280: a = tobf(ARG0(arg),prec);
281: NEWBF(r);
282: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 283: mpfr_cosh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 284: *rp = r;
285: }
286:
287: void mp_tanh(NODE arg,BF *rp)
288: {
289: Num a;
290: int prec;
291: BF r;
292:
293: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
294: a = tobf(ARG0(arg),prec);
295: NEWBF(r);
296: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 297: mpfr_tanh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 298: *rp = r;
299: }
300:
301: void mp_asinh(NODE arg,BF *rp)
302: {
303: Num a;
304: int prec;
305: BF r;
306:
307: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
308: a = tobf(ARG0(arg),prec);
309: NEWBF(r);
310: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 311: mpfr_asinh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 312: *rp = r;
313: }
314: void mp_acosh(NODE arg,BF *rp)
315: {
316: Num a;
317: int prec;
318: BF r;
319:
320: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
321: a = tobf(ARG0(arg),prec);
322: NEWBF(r);
323: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 324: mpfr_acosh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 325: *rp = r;
326: }
327: void mp_atanh(NODE arg,BF *rp)
328: {
329: Num a;
330: int prec;
331: BF r;
332:
333: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
334: a = tobf(ARG0(arg),prec);
335: NEWBF(r);
336: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 337: mpfr_atanh(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 338: *rp = r;
339: }
340:
341: void mp_exp(NODE arg,BF *rp)
342: {
343: Num a;
344: int prec;
345: BF r;
346:
347: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
348: a = tobf(ARG0(arg),prec);
349: NEWBF(r);
350: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 351: mpfr_exp(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 352: *rp = r;
353: }
354:
355: void mp_log(NODE arg,BF *rp)
356: {
357: Num a;
358: int prec;
359: BF r;
360:
361: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
362: a = tobf(ARG0(arg),prec);
363: NEWBF(r);
364: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 365: mpfr_log(r->body,((BF)a)->body,mpfr_roundmode);
1.1 noro 366: *rp = r;
367: }
368:
369: void mp_pow(NODE arg,BF *rp)
370: {
371: Num a,e;
372: int prec;
373: BF r;
374:
375: prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : 0;
376: a = tobf(ARG0(arg),prec);
377: e = tobf(ARG1(arg),prec);
378: NEWBF(r);
379: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 380: mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
1.1 noro 381: *rp = r;
382: }
1.5 noro 383:
384: void Pmpfr_gamma(NODE arg,BF *rp)
385: {
386: Num a;
387: int prec;
388: BF r;
389:
390: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
391: prec *= 3.32193;
392: a = tobf(ARG0(arg),prec);
393: NEWBF(r);
394: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
395: mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
396: *rp = r;
397: }
1.6 takayama 398:
1.7 ! takayama 399: void Pmpfr_floor(NODE arg,Q *rp)
1.6 takayama 400: {
401: Num a;
402: int prec;
403: BF r;
1.7 ! takayama 404: mpz_t t;
! 405: GZ rz;
1.6 takayama 406:
407: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
408: prec *= 3.32193;
409: a = tobf(ARG0(arg),prec);
410: NEWBF(r);
411: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
412: mpfr_floor(r->body,((BF)a)->body);
1.7 ! takayama 413: mpz_init(t);
! 414: mpfr_get_z(t,r->body,mpfr_roundmode);
! 415: MPZTOGZ(t,rz);
! 416: *rp = gztoz(rz);
! 417: }
! 418:
! 419: void Pmpfr_round(NODE arg,Q *rp)
! 420: {
! 421: Num a;
! 422: int prec;
! 423: BF r;
! 424: mpz_t t;
! 425: GZ rz;
! 426:
! 427: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
! 428: prec *= 3.32193;
! 429: a = tobf(ARG0(arg),prec);
! 430: NEWBF(r);
! 431: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
! 432: mpfr_round(r->body,((BF)a)->body);
! 433: mpz_init(t);
! 434: mpfr_get_z(t,r->body,mpfr_roundmode);
! 435: MPZTOGZ(t,rz);
! 436: *rp = gztoz(rz);
1.6 takayama 437: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>