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