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