Annotation of OpenXM_contrib2/asir2000/engine/bf.c, Revision 1.14
1.2 noro 1: /*
1.14 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.13 2015/09/14 05:47:09 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;
1.14 ! noro 95: else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
! 96: (*addnumt[MAX(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;
1.14 ! noro 118: default:
! 119: goto err;
! 120: break;
1.9 noro 121: }
122: MPFRTOBF(r,d);
123: *c = (Num)d;
1.14 ! noro 124: } else if ( NID(b) == N_B ) {
1.9 noro 125: switch ( NID(a) ) {
126: case N_Q:
127: mpfr_init2(r,BFPREC(b));
128: if ( INT((Q)a) ) {
129: z = ztogz((Q)a);
1.11 noro 130: mpfr_add_z(r,((BF)b)->body,z->body,mpfr_roundmode);
1.9 noro 131: } else {
132: q = qtogq((Q)a);
1.11 noro 133: mpfr_add_q(r,((BF)b)->body,q->body,mpfr_roundmode);
1.8 noro 134: }
1.9 noro 135: break;
136: case N_R:
137: /* double precision = 53 */
138: mpfr_init2(r,MAX(BFPREC(b),53));
1.11 noro 139: mpfr_add_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
1.9 noro 140: break;
1.14 ! noro 141: default:
! 142: goto err;
! 143: break;
1.8 noro 144: }
1.7 noro 145: MPFRTOBF(r,d);
146: *c = (Num)d;
1.14 ! noro 147: } else
! 148: goto err;
1.12 noro 149: if ( !cmpbf(*c,0) ) *c = 0;
1.14 ! noro 150: return;
! 151:
! 152: err: error("addbf : invalid argument");
1.7 noro 153: }
154:
155: void subbf(Num a,Num b,Num *c)
156: {
1.8 noro 157: mpfr_t r,s;
158: GZ z;
159: GQ q;
1.7 noro 160: BF d;
161:
162: if ( !a )
163: (*chsgnnumt[NID(b)])(b,c);
164: else if ( !b )
165: *c = a;
1.14 ! noro 166: else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
! 167: (*subnumt[MAX(NID(a),NID(b))])(a,b,c);
1.9 noro 168: else if ( NID(a) == N_B ) {
169: switch ( NID(b) ) {
170: case N_Q:
171: mpfr_init2(r,BFPREC(a));
172: if ( INT((Q)b) ) {
173: z = ztogz((Q)b);
1.11 noro 174: mpfr_sub_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9 noro 175: } else {
176: q = qtogq((Q)b);
1.11 noro 177: mpfr_sub_q(r,((BF)a)->body,q->body,mpfr_roundmode);
1.9 noro 178: }
179: break;
180: case N_R:
181: /* double precision = 53 */
182: mpfr_init2(r,MAX(BFPREC(a),53));
1.11 noro 183: mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
1.9 noro 184: break;
185: case N_B:
186: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11 noro 187: mpfr_sub(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 188: break;
1.14 ! noro 189: default:
! 190: goto err;
1.9 noro 191: }
192: MPFRTOBF(r,d);
193: *c = (Num)d;
1.14 ! noro 194: } else if ( NID(b)==N_B ) {
1.9 noro 195: switch ( NID(a) ) {
196: case N_Q:
197: mpfr_init2(r,BFPREC(b));
198: if ( INT((Q)a) ) {
199: z = ztogz((Q)a);
1.11 noro 200: mpfr_sub_z(r,((BF)b)->body,z->body,mpfr_roundmode);
1.9 noro 201: } else {
202: q = qtogq((Q)a);
1.11 noro 203: mpfr_sub_q(r,((BF)b)->body,q->body,mpfr_roundmode);
1.8 noro 204: }
1.11 noro 205: mpfr_neg(r,r,mpfr_roundmode);
1.9 noro 206: break;
207: case N_R:
208: /* double precision = 53 */
209: mpfr_init2(r,MAX(BFPREC(b),53));
1.11 noro 210: mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 211: break;
1.14 ! noro 212: default:
! 213: goto err;
1.8 noro 214: }
1.14 ! noro 215:
1.7 noro 216: MPFRTOBF(r,d);
217: *c = (Num)d;
1.14 ! noro 218: } else
! 219: goto err;
1.12 noro 220: if ( !cmpbf(*c,0) ) *c = 0;
1.14 ! noro 221: return;
! 222:
! 223: err: error("subbf : invalid argument");
1.7 noro 224: }
225:
226: void mulbf(Num a,Num b,Num *c)
227: {
228: mpfr_t r;
1.8 noro 229: GZ z;
230: GQ q;
1.7 noro 231: BF d;
232: int prec;
233:
234: if ( !a || !b )
235: *c = 0;
1.14 ! noro 236: else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
! 237: (*mulnumt[MAX(NID(a),NID(b))])(a,b,c);
1.9 noro 238: else if ( NID(a) == N_B ) {
239: switch ( NID(b) ) {
240: case N_Q:
241: mpfr_init2(r,BFPREC(a));
242: if ( INT((Q)b) ) {
243: z = ztogz((Q)b);
1.11 noro 244: mpfr_mul_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9 noro 245: } else {
246: q = qtogq((Q)b);
1.11 noro 247: mpfr_mul_q(r,((BF)a)->body,q->body,mpfr_roundmode);
1.9 noro 248: }
249: break;
250: case N_R:
251: /* double precision = 53 */
252: mpfr_init2(r,MAX(BFPREC(a),53));
1.11 noro 253: mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
1.9 noro 254: break;
255: case N_B:
256: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11 noro 257: mpfr_mul(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 258: break;
1.14 ! noro 259: default:
! 260: goto err;
1.9 noro 261: }
262: MPFRTOBF(r,d);
263: *c = (Num)d;
1.14 ! noro 264: } else if ( NID(b) == N_B ) {
1.9 noro 265: switch ( NID(a) ) {
266: case N_Q:
267: mpfr_init2(r,BFPREC(b));
268: if ( INT((Q)a) ) {
269: z = ztogz((Q)a);
1.11 noro 270: mpfr_mul_z(r,((BF)b)->body,z->body,mpfr_roundmode);
1.9 noro 271: } else {
272: q = qtogq((Q)a);
1.11 noro 273: mpfr_mul_q(r,((BF)b)->body,q->body,mpfr_roundmode);
1.8 noro 274: }
1.9 noro 275: break;
276: case N_R:
277: /* double precision = 53 */
278: mpfr_init2(r,MAX(BFPREC(b),53));
1.11 noro 279: mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
1.9 noro 280: break;
1.14 ! noro 281: default:
! 282: goto err;
1.8 noro 283: }
1.7 noro 284: MPFRTOBF(r,d);
285: *c = (Num)d;
1.14 ! noro 286: } else
! 287: goto err;
! 288:
1.12 noro 289: if ( !cmpbf(*c,0) ) *c = 0;
1.14 ! noro 290: return;
! 291:
! 292: err: error("mulbf : invalid argument");
1.7 noro 293: }
294:
295: void divbf(Num a,Num b,Num *c)
296: {
1.8 noro 297: mpfr_t s,r;
298: GZ z;
299: GQ q;
1.7 noro 300: BF d;
301:
302: if ( !b )
303: error("divbf : division by 0");
304: else if ( !a )
305: *c = 0;
1.14 ! noro 306: else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
! 307: (*divnumt[MAX(NID(a),NID(b))])(a,b,c);
1.9 noro 308: else if ( NID(a) == N_B ) {
309: switch ( NID(b) ) {
310: case N_Q:
311: mpfr_init2(r,BFPREC(a));
312: if ( INT((Q)b) ) {
313: z = ztogz((Q)b);
1.11 noro 314: mpfr_div_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9 noro 315: } else {
316: q = qtogq((Q)b);
1.11 noro 317: mpfr_div_q(r,((BF)a)->body,q->body,mpfr_roundmode);
1.8 noro 318: }
1.9 noro 319: break;
320: case N_R:
321: /* double precision = 53 */
322: mpfr_init2(r,MAX(BFPREC(a),53));
1.11 noro 323: mpfr_div_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
1.9 noro 324: break;
325: case N_B:
326: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11 noro 327: mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 328: break;
1.14 ! noro 329: default:
! 330: goto err;
1.9 noro 331: }
332: MPFRTOBF(r,d);
333: *c = (Num)d;
1.14 ! noro 334: } else if ( NID(b)==N_B ) {
1.9 noro 335: switch ( NID(a) ) {
336: case N_Q:
337: /* XXX : mpfr_z_div and mpfr_q_div are not implemented */
338: a = tobf(a,BFPREC(b));
339: mpfr_init2(r,BFPREC(b));
1.11 noro 340: mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 341: break;
342: case N_R:
343: /* double precision = 53 */
344: mpfr_init2(r,MAX(BFPREC(b),53));
1.11 noro 345: mpfr_d_div(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 346: break;
1.14 ! noro 347: default:
! 348: goto err;
1.8 noro 349: }
1.7 noro 350: MPFRTOBF(r,d);
351: *c = (Num)d;
1.14 ! noro 352: } else
! 353: goto err;
! 354:
1.12 noro 355: if ( !cmpbf(*c,0) ) *c = 0;
1.14 ! noro 356: return;
! 357:
! 358: err: error("mulbf : invalid argument");
1.7 noro 359: }
360:
361: void pwrbf(Num a,Num b,Num *c)
362: {
1.8 noro 363: int prec;
1.7 noro 364: mpfr_t r;
1.8 noro 365: GZ z;
1.7 noro 366: BF d;
367:
368: if ( !b )
369: *c = (Num)ONE;
370: else if ( !a )
371: *c = 0;
1.14 ! noro 372: else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
! 373: (*pwrnumt[MAX(NID(a),NID(b))])(a,b,c);
1.9 noro 374: else if ( NID(a) == N_B ) {
375: switch ( NID(b) ) {
376: case N_Q:
377: mpfr_init2(r,BFPREC(a));
378: if ( INT((Q)b) ) {
379: z = ztogz((Q)b);
1.11 noro 380: mpfr_pow_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9 noro 381: } else {
382: b = tobf(b,BFPREC(a));
1.11 noro 383: mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.8 noro 384: }
1.9 noro 385: break;
386: case N_R:
387: /* double precision = 53 */
388: prec = MAX(BFPREC(a),53);
389: mpfr_init2(r,prec);
390: b = tobf(b,prec);
1.11 noro 391: mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 392: break;
393: case N_B:
394: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11 noro 395: mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 396: break;
1.14 ! noro 397: default:
! 398: goto err;
1.9 noro 399: }
400: MPFRTOBF(r,d);
401: *c = (Num)d;
1.14 ! noro 402: } else if ( NID(b)==N_B ) {
1.9 noro 403: switch ( NID(a) ) {
404: case N_Q:
405: mpfr_init2(r,BFPREC(b));
406: a = tobf(a,BFPREC(b));
1.11 noro 407: mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 408: break;
409: case N_R:
410: /* double precision = 53 */
411: prec = MAX(BFPREC(a),53);
412: mpfr_init2(r,prec);
413: a = tobf(a,prec);
1.11 noro 414: mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9 noro 415: break;
1.14 ! noro 416: default:
! 417: goto err;
1.8 noro 418: }
1.7 noro 419: MPFRTOBF(r,d);
420: *c = (Num)d;
1.14 ! noro 421: } else
! 422: goto err;
! 423:
1.12 noro 424: if ( !cmpbf(*c,0) ) *c = 0;
1.14 ! noro 425: return;
! 426:
! 427: err: error("pwrbf : invalid argument");
1.7 noro 428: }
429:
430: void chsgnbf(Num a,Num *c)
431: {
432: mpfr_t r;
433: BF d;
434:
435: if ( !a )
436: *c = 0;
1.14 ! noro 437: else if ( NID(a) <= N_R )
1.7 noro 438: (*chsgnnumt[NID(a)])(a,c);
1.14 ! noro 439: else if ( NID(a) == N_B ) {
1.7 noro 440: mpfr_init2(r,BFPREC(a));
1.11 noro 441: mpfr_neg(r,((BF)a)->body,mpfr_roundmode);
1.7 noro 442: MPFRTOBF(r,d);
443: *c = (Num)d;
1.14 ! noro 444: } else
! 445: error("chsgnbf : invalid argument");
1.7 noro 446: }
447:
448: int cmpbf(Num a,Num b)
449: {
1.9 noro 450: int ret;
451: GZ z;
452: GQ q;
453:
1.7 noro 454: if ( !a ) {
1.12 noro 455: if ( !b ) return 0;
1.14 ! noro 456: else if ( NID(b)<=N_R )
1.7 noro 457: return (*cmpnumt[NID(b)])(a,b);
1.14 ! noro 458: else if ( NID(b)==N_B )
! 459: return -mpfr_sgn(((BF)b)->body);
1.7 noro 460: else
1.14 ! noro 461: goto err;
1.7 noro 462: } else if ( !b ) {
1.14 ! noro 463: if ( NID(a)<=N_R )
1.7 noro 464: return (*cmpnumt[NID(a)])(a,b);
1.14 ! noro 465: else if ( NID(a)==N_B )
! 466: return mpfr_sgn(((BF)a)->body);
1.7 noro 467: else
1.14 ! noro 468: goto err;
! 469: } else if ( NID(a) <= N_R && NID(b) <= N_R )
! 470: return (*cmpnumt[MAX(NID(a),NID(b))])(a,b);
! 471: else if ( NID(a) == N_B ) {
1.9 noro 472: switch ( NID(b) ) {
473: case N_Q:
474: if ( INT((Q)b) ) {
475: z = ztogz((Q)b);
476: ret = mpfr_cmp_z(((BF)a)->body,z->body);
477: } else {
478: q = qtogq((Q)b);
479: ret = mpfr_cmp_q(((BF)a)->body,q->body);
480: }
481: break;
482: case N_R:
483: /* double precision = 53 */
484: ret = mpfr_cmp_d(((BF)a)->body,((Real)b)->body);
485: break;
486: case N_B:
487: ret = mpfr_cmp(((BF)a)->body,((BF)b)->body);
488: break;
1.14 ! noro 489: default:
! 490: goto err;
1.9 noro 491: }
492: return ret;
1.14 ! noro 493: } else if ( NID(b)==N_B ) {
1.9 noro 494: switch ( NID(a) ) {
495: case N_Q:
496: if ( INT((Q)a) ) {
497: z = ztogz((Q)a);
498: ret = mpfr_cmp_z(((BF)b)->body,z->body);
499: } else {
500: q = qtogq((Q)a);
501: ret = mpfr_cmp_q(((BF)b)->body,q->body);
502: }
503: break;
504: case N_R:
505: /* double precision = 53 */
506: ret = mpfr_cmp_d(((BF)b)->body,((Real)a)->body);
507: break;
1.14 ! noro 508: default:
! 509: goto err;
1.9 noro 510: }
511: return -ret;
1.7 noro 512: }
1.14 ! noro 513: err: error("cmpbf : cannot compare");
1.1 noro 514: }
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