Annotation of OpenXM_contrib2/asir2000/engine/bf.c, Revision 1.9
1.2 noro 1: /*
1.9 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.8 2015/08/05 03:46:35 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(s,BFPREC(b));
! 176: mpfr_init2(r,BFPREC(b));
! 177: if ( INT((Q)a) ) {
! 178: z = ztogz((Q)a);
! 179: mpfr_sub_z(s,((BF)b)->body,z->body,MPFR_RNDN);
! 180: } else {
! 181: q = qtogq((Q)a);
! 182: mpfr_sub_q(s,((BF)b)->body,q->body,MPFR_RNDN);
1.8 noro 183: }
1.9 ! noro 184: mpfr_neg(r,s,MPFR_RNDN);
! 185: break;
! 186: case N_R:
! 187: /* double precision = 53 */
! 188: mpfr_init2(r,MAX(BFPREC(b),53));
! 189: mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);
! 190: break;
1.8 noro 191: }
1.7 noro 192: MPFRTOBF(r,d);
193: *c = (Num)d;
194: }
195: }
196:
197: void mulbf(Num a,Num b,Num *c)
198: {
199: mpfr_t r;
1.8 noro 200: GZ z;
201: GQ q;
1.7 noro 202: BF d;
203: int prec;
204:
205: if ( !a || !b )
206: *c = 0;
207: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
208: (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9 ! noro 209: else if ( NID(a) == N_B ) {
! 210: switch ( NID(b) ) {
! 211: case N_Q:
! 212: mpfr_init2(r,BFPREC(a));
! 213: if ( INT((Q)b) ) {
! 214: z = ztogz((Q)b);
! 215: mpfr_mul_z(r,((BF)a)->body,z->body,MPFR_RNDN);
! 216: } else {
! 217: q = qtogq((Q)b);
! 218: mpfr_mul_q(r,((BF)a)->body,q->body,MPFR_RNDN);
! 219: }
! 220: break;
! 221: case N_R:
! 222: /* double precision = 53 */
! 223: mpfr_init2(r,MAX(BFPREC(a),53));
! 224: mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
! 225: break;
! 226: case N_B:
! 227: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
! 228: mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 229: break;
! 230: }
! 231: MPFRTOBF(r,d);
! 232: *c = (Num)d;
! 233: } else { /* NID(b)==N_B */
! 234: switch ( NID(a) ) {
! 235: case N_Q:
! 236: mpfr_init2(r,BFPREC(b));
! 237: if ( INT((Q)a) ) {
! 238: z = ztogz((Q)a);
! 239: mpfr_mul_z(r,((BF)b)->body,z->body,MPFR_RNDN);
! 240: } else {
! 241: q = qtogq((Q)a);
! 242: mpfr_mul_q(r,((BF)b)->body,q->body,MPFR_RNDN);
1.8 noro 243: }
1.9 ! noro 244: break;
! 245: case N_R:
! 246: /* double precision = 53 */
! 247: mpfr_init2(r,MAX(BFPREC(b),53));
! 248: mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);
! 249: break;
1.8 noro 250: }
1.7 noro 251: MPFRTOBF(r,d);
252: *c = (Num)d;
253: }
254: }
255:
256: void divbf(Num a,Num b,Num *c)
257: {
1.8 noro 258: mpfr_t s,r;
259: GZ z;
260: GQ q;
1.7 noro 261: BF d;
262:
263: if ( !b )
264: error("divbf : division by 0");
265: else if ( !a )
266: *c = 0;
267: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
268: (*divnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9 ! noro 269: else if ( NID(a) == N_B ) {
! 270: switch ( NID(b) ) {
! 271: case N_Q:
! 272: mpfr_init2(r,BFPREC(a));
! 273: if ( INT((Q)b) ) {
! 274: z = ztogz((Q)b);
! 275: mpfr_div_z(r,((BF)a)->body,z->body,MPFR_RNDN);
! 276: } else {
! 277: q = qtogq((Q)b);
! 278: mpfr_div_q(r,((BF)a)->body,q->body,MPFR_RNDN);
1.8 noro 279: }
1.9 ! noro 280: break;
! 281: case N_R:
! 282: /* double precision = 53 */
! 283: mpfr_init2(r,MAX(BFPREC(a),53));
! 284: mpfr_div_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
! 285: break;
! 286: case N_B:
! 287: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
! 288: mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 289: break;
! 290: }
! 291: MPFRTOBF(r,d);
! 292: *c = (Num)d;
! 293: } else { /* NID(b)==N_B */
! 294: switch ( NID(a) ) {
! 295: case N_Q:
! 296: /* XXX : mpfr_z_div and mpfr_q_div are not implemented */
! 297: a = tobf(a,BFPREC(b));
! 298: mpfr_init2(r,BFPREC(b));
! 299: mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 300: break;
! 301: case N_R:
! 302: /* double precision = 53 */
! 303: mpfr_init2(r,MAX(BFPREC(b),53));
! 304: mpfr_d_div(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);
! 305: break;
1.8 noro 306: }
1.7 noro 307: MPFRTOBF(r,d);
308: *c = (Num)d;
309: }
310: }
311:
312: void pwrbf(Num a,Num b,Num *c)
313: {
1.8 noro 314: int prec;
1.7 noro 315: mpfr_t r;
1.8 noro 316: GZ z;
1.7 noro 317: BF d;
318:
319: if ( !b )
320: *c = (Num)ONE;
321: else if ( !a )
322: *c = 0;
323: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
324: (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9 ! noro 325: else if ( NID(a) == N_B ) {
! 326: switch ( NID(b) ) {
! 327: case N_Q:
! 328: mpfr_init2(r,BFPREC(a));
! 329: if ( INT((Q)b) ) {
! 330: z = ztogz((Q)b);
! 331: mpfr_pow_z(r,((BF)a)->body,z->body,MPFR_RNDN);
! 332: } else {
! 333: b = tobf(b,BFPREC(a));
1.8 noro 334: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
335: }
1.9 ! noro 336: break;
! 337: case N_R:
! 338: /* double precision = 53 */
! 339: prec = MAX(BFPREC(a),53);
! 340: mpfr_init2(r,prec);
! 341: b = tobf(b,prec);
! 342: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 343: break;
! 344: case N_B:
! 345: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
! 346: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 347: break;
! 348: }
! 349: MPFRTOBF(r,d);
! 350: *c = (Num)d;
! 351: } else { /* NID(b)==N_B */
! 352: switch ( NID(a) ) {
! 353: case N_Q:
! 354: mpfr_init2(r,BFPREC(b));
! 355: a = tobf(a,BFPREC(b));
! 356: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 357: break;
! 358: case N_R:
! 359: /* double precision = 53 */
! 360: prec = MAX(BFPREC(a),53);
! 361: mpfr_init2(r,prec);
! 362: a = tobf(a,prec);
! 363: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 364: break;
1.8 noro 365: }
1.7 noro 366: MPFRTOBF(r,d);
367: *c = (Num)d;
368: }
369: }
370:
371: void chsgnbf(Num a,Num *c)
372: {
373: mpfr_t r;
374: BF d;
375:
376: if ( !a )
377: *c = 0;
378: else if ( NID(a) <= N_A )
379: (*chsgnnumt[NID(a)])(a,c);
380: else {
381: mpfr_init2(r,BFPREC(a));
382: mpfr_neg(r,((BF)a)->body,MPFR_RNDN);
383: MPFRTOBF(r,d);
384: *c = (Num)d;
385: }
386: }
387:
388: int cmpbf(Num a,Num b)
389: {
1.9 ! noro 390: int ret;
! 391: GZ z;
! 392: GQ q;
! 393:
1.7 noro 394: if ( !a ) {
395: if ( !b || (NID(b)<=N_A) )
396: return (*cmpnumt[NID(b)])(a,b);
397: else
398: return -mpfr_sgn(((BF)a)->body);
399: } else if ( !b ) {
400: if ( !a || (NID(a)<=N_A) )
401: return (*cmpnumt[NID(a)])(a,b);
402: else
403: return mpfr_sgn(((BF)a)->body);
1.9 ! noro 404: } else if ( NID(a) == N_B ) {
! 405: switch ( NID(b) ) {
! 406: case N_Q:
! 407: if ( INT((Q)b) ) {
! 408: z = ztogz((Q)b);
! 409: ret = mpfr_cmp_z(((BF)a)->body,z->body);
! 410: } else {
! 411: q = qtogq((Q)b);
! 412: ret = mpfr_cmp_q(((BF)a)->body,q->body);
! 413: }
! 414: break;
! 415: case N_R:
! 416: /* double precision = 53 */
! 417: ret = mpfr_cmp_d(((BF)a)->body,((Real)b)->body);
! 418: break;
! 419: case N_B:
! 420: ret = mpfr_cmp(((BF)a)->body,((BF)b)->body);
! 421: break;
! 422: }
! 423: return ret;
! 424: } else { /* NID(b)==N_B */
! 425: switch ( NID(a) ) {
! 426: case N_Q:
! 427: if ( INT((Q)a) ) {
! 428: z = ztogz((Q)a);
! 429: ret = mpfr_cmp_z(((BF)b)->body,z->body);
! 430: } else {
! 431: q = qtogq((Q)a);
! 432: ret = mpfr_cmp_q(((BF)b)->body,q->body);
! 433: }
! 434: break;
! 435: case N_R:
! 436: /* double precision = 53 */
! 437: ret = mpfr_cmp_d(((BF)b)->body,((Real)a)->body);
! 438: break;
! 439: }
! 440: return -ret;
1.7 noro 441: }
1.1 noro 442: }
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