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