Annotation of OpenXM_contrib2/asir2000/engine/bf.c, Revision 1.10
1.2 noro 1: /*
1.10 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.9 2015/08/06 00:43:20 noro Exp $
1.7 noro 3: */
1.1 noro 4: #include "ca.h"
5: #include "base.h"
6: #include <math.h>
7:
1.7 noro 8: Num tobf(Num,int);
1.1 noro 9:
1.7 noro 10: #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)
11:
12: void strtobf(char *s,BF *p)
13: {
14: BF r;
15: NEWBF(r);
16: mpfr_init(r->body);
17: mpfr_set_str(r->body,s,10,MPFR_RNDN);
18: *p = r;
19: }
20:
21: double mpfrtodbl(mpfr_t a)
22: {
23: return mpfr_get_d(a,MPFR_RNDN);
24: }
25:
26: Num tobf(Num a,int prec)
27: {
28: mpfr_t r;
29: mpz_t z;
30: mpq_t q;
31: BF d;
32: N nm,dn;
33: int sgn;
34:
35: if ( !a ) {
1.8 noro 36: prec ? mpfr_init2(r,prec) : mpfr_init(r);
1.7 noro 37: mpfr_set_zero(r,1);
38: MPFRTOBF(r,d);
39: return (Num)d;
40: } else {
41: switch ( NID(a) ) {
42: case N_B:
43: return a;
44: break;
1.8 noro 45: case N_R:
46: prec ? mpfr_init2(r,prec) : mpfr_init(r);
47: mpfr_init_set_d(r,((Real)a)->body,MPFR_RNDN);
1.7 noro 48: MPFRTOBF(r,d);
49: return (Num)d;
1.8 noro 50: break;
1.7 noro 51: case N_Q:
52: nm = NM((Q)a); dn = DN((Q)a); sgn = SGN((Q)a);
53: if ( INT((Q)a) ) {
54: mpz_init(z);
55: mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
56: if ( sgn < 0 ) mpz_neg(z,z);
57: mpfr_init_set_z(r,z,MPFR_RNDN);
58: } else {
59: mpq_init(q);
60: mpz_import(mpq_numref(q),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
61: mpz_import(mpq_denref(q),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));
62: if ( sgn < 0 ) mpq_neg(q,q);
63: mpfr_init_set_q(r,q,MPFR_RNDN);
64: }
65: MPFRTOBF(r,d);
66: return (Num)d;
67: break;
68: default:
69: error("tobf : invalid argument");
1.8 noro 70: break;
1.7 noro 71: }
72: }
73: }
74:
75: void addbf(Num a,Num b,Num *c)
76: {
77: mpfr_t r;
78: BF d;
1.8 noro 79: GZ z;
80: GQ q;
1.7 noro 81:
82: if ( !a )
83: *c = b;
84: else if ( !b )
85: *c = a;
86: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
87: (*addnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9 noro 88: else if ( NID(a) == N_B ) {
89: switch ( NID(b) ) {
90: case N_Q:
91: mpfr_init2(r,BFPREC(a));
92: if ( INT((Q)b) ) {
93: z = ztogz((Q)b);
94: mpfr_add_z(r,((BF)a)->body,z->body,MPFR_RNDN);
95: } else {
96: q = qtogq((Q)b);
97: mpfr_add_q(r,((BF)a)->body,q->body,MPFR_RNDN);
98: }
99: break;
100: case N_R:
101: /* double precision = 53 */
102: mpfr_init2(r,MAX(BFPREC(a),53));
103: mpfr_add_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
104: break;
105: case N_B:
106: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
107: mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
108: break;
109: }
110: MPFRTOBF(r,d);
111: *c = (Num)d;
112: } else { /* NID(b)==N_B */
113: switch ( NID(a) ) {
114: case N_Q:
115: mpfr_init2(r,BFPREC(b));
116: if ( INT((Q)a) ) {
117: z = ztogz((Q)a);
118: mpfr_add_z(r,((BF)b)->body,z->body,MPFR_RNDN);
119: } else {
120: q = qtogq((Q)a);
121: mpfr_add_q(r,((BF)b)->body,q->body,MPFR_RNDN);
1.8 noro 122: }
1.9 noro 123: break;
124: case N_R:
125: /* double precision = 53 */
126: mpfr_init2(r,MAX(BFPREC(b),53));
127: mpfr_add_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);
128: break;
1.8 noro 129: }
1.7 noro 130: MPFRTOBF(r,d);
131: *c = (Num)d;
132: }
133: }
134:
135: void subbf(Num a,Num b,Num *c)
136: {
1.8 noro 137: mpfr_t r,s;
138: GZ z;
139: GQ q;
1.7 noro 140: BF d;
141:
142: if ( !a )
143: (*chsgnnumt[NID(b)])(b,c);
144: else if ( !b )
145: *c = a;
146: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
147: (*subnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9 noro 148: else if ( NID(a) == N_B ) {
149: switch ( NID(b) ) {
150: case N_Q:
151: mpfr_init2(r,BFPREC(a));
152: if ( INT((Q)b) ) {
153: z = ztogz((Q)b);
154: mpfr_sub_z(r,((BF)a)->body,z->body,MPFR_RNDN);
155: } else {
156: q = qtogq((Q)b);
157: mpfr_sub_q(r,((BF)a)->body,q->body,MPFR_RNDN);
158: }
159: break;
160: case N_R:
161: /* double precision = 53 */
162: mpfr_init2(r,MAX(BFPREC(a),53));
163: mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
164: break;
165: case N_B:
166: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
167: mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
168: break;
169: }
170: MPFRTOBF(r,d);
171: *c = (Num)d;
172: } else { /* NID(b)==N_B */
173: switch ( NID(a) ) {
174: case N_Q:
175: mpfr_init2(r,BFPREC(b));
176: if ( INT((Q)a) ) {
177: z = ztogz((Q)a);
1.10 ! noro 178: mpfr_sub_z(r,((BF)b)->body,z->body,MPFR_RNDN);
1.9 noro 179: } else {
180: q = qtogq((Q)a);
1.10 ! noro 181: mpfr_sub_q(r,((BF)b)->body,q->body,MPFR_RNDN);
1.8 noro 182: }
1.10 ! noro 183: mpfr_neg(r,r,MPFR_RNDN);
1.9 noro 184: break;
185: case N_R:
186: /* double precision = 53 */
187: mpfr_init2(r,MAX(BFPREC(b),53));
188: mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);
189: break;
1.8 noro 190: }
1.7 noro 191: MPFRTOBF(r,d);
192: *c = (Num)d;
193: }
194: }
195:
196: void mulbf(Num a,Num b,Num *c)
197: {
198: mpfr_t r;
1.8 noro 199: GZ z;
200: GQ q;
1.7 noro 201: BF d;
202: int prec;
203:
204: if ( !a || !b )
205: *c = 0;
206: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
207: (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9 noro 208: else if ( NID(a) == N_B ) {
209: switch ( NID(b) ) {
210: case N_Q:
211: mpfr_init2(r,BFPREC(a));
212: if ( INT((Q)b) ) {
213: z = ztogz((Q)b);
214: mpfr_mul_z(r,((BF)a)->body,z->body,MPFR_RNDN);
215: } else {
216: q = qtogq((Q)b);
217: mpfr_mul_q(r,((BF)a)->body,q->body,MPFR_RNDN);
218: }
219: break;
220: case N_R:
221: /* double precision = 53 */
222: mpfr_init2(r,MAX(BFPREC(a),53));
223: mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
224: break;
225: case N_B:
226: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
227: mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
228: break;
229: }
230: MPFRTOBF(r,d);
231: *c = (Num)d;
232: } else { /* NID(b)==N_B */
233: switch ( NID(a) ) {
234: case N_Q:
235: mpfr_init2(r,BFPREC(b));
236: if ( INT((Q)a) ) {
237: z = ztogz((Q)a);
238: mpfr_mul_z(r,((BF)b)->body,z->body,MPFR_RNDN);
239: } else {
240: q = qtogq((Q)a);
241: mpfr_mul_q(r,((BF)b)->body,q->body,MPFR_RNDN);
1.8 noro 242: }
1.9 noro 243: break;
244: case N_R:
245: /* double precision = 53 */
246: mpfr_init2(r,MAX(BFPREC(b),53));
247: mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);
248: break;
1.8 noro 249: }
1.7 noro 250: MPFRTOBF(r,d);
251: *c = (Num)d;
252: }
253: }
254:
255: void divbf(Num a,Num b,Num *c)
256: {
1.8 noro 257: mpfr_t s,r;
258: GZ z;
259: GQ q;
1.7 noro 260: BF d;
261:
262: if ( !b )
263: error("divbf : division by 0");
264: else if ( !a )
265: *c = 0;
266: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
267: (*divnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9 noro 268: else if ( NID(a) == N_B ) {
269: switch ( NID(b) ) {
270: case N_Q:
271: mpfr_init2(r,BFPREC(a));
272: if ( INT((Q)b) ) {
273: z = ztogz((Q)b);
274: mpfr_div_z(r,((BF)a)->body,z->body,MPFR_RNDN);
275: } else {
276: q = qtogq((Q)b);
277: mpfr_div_q(r,((BF)a)->body,q->body,MPFR_RNDN);
1.8 noro 278: }
1.9 noro 279: break;
280: case N_R:
281: /* double precision = 53 */
282: mpfr_init2(r,MAX(BFPREC(a),53));
283: mpfr_div_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
284: break;
285: case N_B:
286: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
287: mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
288: break;
289: }
290: MPFRTOBF(r,d);
291: *c = (Num)d;
292: } else { /* NID(b)==N_B */
293: switch ( NID(a) ) {
294: case N_Q:
295: /* XXX : mpfr_z_div and mpfr_q_div are not implemented */
296: a = tobf(a,BFPREC(b));
297: mpfr_init2(r,BFPREC(b));
298: mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
299: break;
300: case N_R:
301: /* double precision = 53 */
302: mpfr_init2(r,MAX(BFPREC(b),53));
303: mpfr_d_div(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);
304: break;
1.8 noro 305: }
1.7 noro 306: MPFRTOBF(r,d);
307: *c = (Num)d;
308: }
309: }
310:
311: void pwrbf(Num a,Num b,Num *c)
312: {
1.8 noro 313: int prec;
1.7 noro 314: mpfr_t r;
1.8 noro 315: GZ z;
1.7 noro 316: BF d;
317:
318: if ( !b )
319: *c = (Num)ONE;
320: else if ( !a )
321: *c = 0;
322: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
323: (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9 noro 324: else if ( NID(a) == N_B ) {
325: switch ( NID(b) ) {
326: case N_Q:
327: mpfr_init2(r,BFPREC(a));
328: if ( INT((Q)b) ) {
329: z = ztogz((Q)b);
330: mpfr_pow_z(r,((BF)a)->body,z->body,MPFR_RNDN);
331: } else {
332: b = tobf(b,BFPREC(a));
1.8 noro 333: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
334: }
1.9 noro 335: break;
336: case N_R:
337: /* double precision = 53 */
338: prec = MAX(BFPREC(a),53);
339: mpfr_init2(r,prec);
340: b = tobf(b,prec);
341: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
342: break;
343: case N_B:
344: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
345: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
346: break;
347: }
348: MPFRTOBF(r,d);
349: *c = (Num)d;
350: } else { /* NID(b)==N_B */
351: switch ( NID(a) ) {
352: case N_Q:
353: mpfr_init2(r,BFPREC(b));
354: a = tobf(a,BFPREC(b));
355: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
356: break;
357: case N_R:
358: /* double precision = 53 */
359: prec = MAX(BFPREC(a),53);
360: mpfr_init2(r,prec);
361: a = tobf(a,prec);
362: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
363: break;
1.8 noro 364: }
1.7 noro 365: MPFRTOBF(r,d);
366: *c = (Num)d;
367: }
368: }
369:
370: void chsgnbf(Num a,Num *c)
371: {
372: mpfr_t r;
373: BF d;
374:
375: if ( !a )
376: *c = 0;
377: else if ( NID(a) <= N_A )
378: (*chsgnnumt[NID(a)])(a,c);
379: else {
380: mpfr_init2(r,BFPREC(a));
381: mpfr_neg(r,((BF)a)->body,MPFR_RNDN);
382: MPFRTOBF(r,d);
383: *c = (Num)d;
384: }
385: }
386:
387: int cmpbf(Num a,Num b)
388: {
1.9 noro 389: int ret;
390: GZ z;
391: GQ q;
392:
1.7 noro 393: if ( !a ) {
394: if ( !b || (NID(b)<=N_A) )
395: return (*cmpnumt[NID(b)])(a,b);
396: else
397: return -mpfr_sgn(((BF)a)->body);
398: } else if ( !b ) {
399: if ( !a || (NID(a)<=N_A) )
400: return (*cmpnumt[NID(a)])(a,b);
401: else
402: return mpfr_sgn(((BF)a)->body);
1.9 noro 403: } else if ( NID(a) == N_B ) {
404: switch ( NID(b) ) {
405: case N_Q:
406: if ( INT((Q)b) ) {
407: z = ztogz((Q)b);
408: ret = mpfr_cmp_z(((BF)a)->body,z->body);
409: } else {
410: q = qtogq((Q)b);
411: ret = mpfr_cmp_q(((BF)a)->body,q->body);
412: }
413: break;
414: case N_R:
415: /* double precision = 53 */
416: ret = mpfr_cmp_d(((BF)a)->body,((Real)b)->body);
417: break;
418: case N_B:
419: ret = mpfr_cmp(((BF)a)->body,((BF)b)->body);
420: break;
421: }
422: return ret;
423: } else { /* NID(b)==N_B */
424: switch ( NID(a) ) {
425: case N_Q:
426: if ( INT((Q)a) ) {
427: z = ztogz((Q)a);
428: ret = mpfr_cmp_z(((BF)b)->body,z->body);
429: } else {
430: q = qtogq((Q)a);
431: ret = mpfr_cmp_q(((BF)b)->body,q->body);
432: }
433: break;
434: case N_R:
435: /* double precision = 53 */
436: ret = mpfr_cmp_d(((BF)b)->body,((Real)a)->body);
437: break;
438: }
439: return -ret;
1.7 noro 440: }
1.1 noro 441: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bf.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine