Annotation of OpenXM_contrib2/asir2000/engine/bf.c, Revision 1.8
1.2 noro 1: /*
1.8 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.7 2015/08/04 06:20:45 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);
88: else {
1.8 ! noro 89: if ( NID(a) == N_B ) {
! 90: switch ( NID(b) ) {
! 91: case N_Q:
! 92: mpfr_init2(r,BFPREC(a));
! 93: if ( INT((Q)b) ) {
! 94: z = ztogz((Q)b);
! 95: mpfr_add_z(r,((BF)a)->body,z->body,MPFR_RNDN);
! 96: } else {
! 97: q = qtogq((Q)b);
! 98: mpfr_add_q(r,((BF)a)->body,q->body,MPFR_RNDN);
! 99: }
! 100: break;
! 101: case N_R:
! 102: /* double precision = 53 */
! 103: mpfr_init2(r,MIN(BFPREC(a),53));
! 104: mpfr_add_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
! 105: break;
! 106: case N_B:
! 107: mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
! 108: mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 109: break;
! 110: }
! 111: } else { /* NID(b)==N_B */
! 112: switch ( NID(a) ) {
! 113: case N_Q:
! 114: mpfr_init2(r,BFPREC(b));
! 115: if ( INT((Q)a) ) {
! 116: z = ztogz((Q)a);
! 117: mpfr_add_z(r,((BF)b)->body,z->body,MPFR_RNDN);
! 118: } else {
! 119: q = qtogq((Q)a);
! 120: mpfr_add_q(r,((BF)b)->body,q->body,MPFR_RNDN);
! 121: }
! 122: break;
! 123: case N_R:
! 124: /* double precision = 53 */
! 125: mpfr_init2(r,MIN(BFPREC(b),53));
! 126: mpfr_add_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);
! 127: break;
! 128: }
! 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);
148: else {
1.8 ! noro 149: if ( NID(a) == N_B ) {
! 150: switch ( NID(b) ) {
! 151: case N_Q:
! 152: mpfr_init2(r,BFPREC(a));
! 153: if ( INT((Q)b) ) {
! 154: z = ztogz((Q)b);
! 155: mpfr_sub_z(r,((BF)a)->body,z->body,MPFR_RNDN);
! 156: } else {
! 157: q = qtogq((Q)b);
! 158: mpfr_sub_q(r,((BF)a)->body,q->body,MPFR_RNDN);
! 159: }
! 160: break;
! 161: case N_R:
! 162: /* double precision = 53 */
! 163: mpfr_init2(r,MIN(BFPREC(a),53));
! 164: mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
! 165: break;
! 166: case N_B:
! 167: mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
! 168: mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 169: break;
! 170: }
! 171: } else { /* NID(b)==N_B */
! 172: switch ( NID(a) ) {
! 173: case N_Q:
! 174: mpfr_init2(s,BFPREC(b));
! 175: mpfr_init2(r,BFPREC(b));
! 176: if ( INT((Q)a) ) {
! 177: z = ztogz((Q)a);
! 178: mpfr_sub_z(s,((BF)b)->body,z->body,MPFR_RNDN);
! 179: } else {
! 180: q = qtogq((Q)a);
! 181: mpfr_sub_q(s,((BF)b)->body,q->body,MPFR_RNDN);
! 182: }
! 183: mpfr_neg(r,s,MPFR_RNDN);
! 184: break;
! 185: case N_R:
! 186: /* double precision = 53 */
! 187: mpfr_init2(r,MIN(BFPREC(b),53));
! 188: mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);
! 189: break;
! 190: }
! 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);
209: else {
1.8 ! noro 210: 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);
! 216: mpfr_mul_z(r,((BF)a)->body,z->body,MPFR_RNDN);
! 217: } else {
! 218: q = qtogq((Q)b);
! 219: mpfr_mul_q(r,((BF)a)->body,q->body,MPFR_RNDN);
! 220: }
! 221: break;
! 222: case N_R:
! 223: /* double precision = 53 */
! 224: mpfr_init2(r,MIN(BFPREC(a),53));
! 225: mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
! 226: break;
! 227: case N_B:
! 228: mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
! 229: mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 230: break;
! 231: }
! 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);
! 242: }
! 243: break;
! 244: case N_R:
! 245: /* double precision = 53 */
! 246: mpfr_init2(r,MIN(BFPREC(b),53));
! 247: mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);
! 248: break;
! 249: }
! 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);
269: else {
1.8 ! noro 270: 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);
! 276: mpfr_div_z(r,((BF)a)->body,z->body,MPFR_RNDN);
! 277: } else {
! 278: q = qtogq((Q)b);
! 279: mpfr_div_q(r,((BF)a)->body,q->body,MPFR_RNDN);
! 280: }
! 281: break;
! 282: case N_R:
! 283: /* double precision = 53 */
! 284: mpfr_init2(r,MIN(BFPREC(a),53));
! 285: mpfr_div_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
! 286: break;
! 287: case N_B:
! 288: mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
! 289: mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 290: break;
! 291: }
! 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,MIN(BFPREC(b),53));
! 303: mpfr_d_div(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);
! 304: break;
! 305: }
! 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);
325: else {
1.8 ! noro 326: 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);
! 332: mpfr_pow_z(r,((BF)a)->body,z->body,MPFR_RNDN);
! 333: } else {
! 334: b = tobf(b,BFPREC(a));
! 335: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 336: }
! 337: break;
! 338: case N_R:
! 339: /* double precision = 53 */
! 340: prec = MIN(BFPREC(a),53);
! 341: mpfr_init2(r,prec);
! 342: b = tobf(b,prec);
! 343: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 344: break;
! 345: case N_B:
! 346: mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
! 347: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 348: break;
! 349: }
! 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 = MIN(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;
! 364: }
! 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: {
390: if ( !a ) {
391: if ( !b || (NID(b)<=N_A) )
392: return (*cmpnumt[NID(b)])(a,b);
393: else
394: return -mpfr_sgn(((BF)a)->body);
395: } else if ( !b ) {
396: if ( !a || (NID(a)<=N_A) )
397: return (*cmpnumt[NID(a)])(a,b);
398: else
399: return mpfr_sgn(((BF)a)->body);
400: } else {
401: if ( NID(a) != N_B ) a = tobf(a,0);
402: if ( NID(b) != N_B ) b = tobf(b,0);
403: return mpfr_cmp(((BF)a)->body,((BF)b)->body);
404: }
1.1 noro 405: }
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