Annotation of OpenXM_contrib2/asir2000/engine/bf.c, Revision 1.7
1.2 noro 1: /*
1.7 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.6 2009/07/15 02:19:47 noro Exp $
! 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 ) {
! 36: prec ? mpfr_init2(r,prec) : mpfr_init(r);
! 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;
! 45: case N_R:
! 46: prec ? mpfr_init2(r,prec) : mpfr_init(r);
! 47: mpfr_init_set_d(r,((Real)a)->body,MPFR_RNDN);
! 48: MPFRTOBF(r,d);
! 49: return (Num)d;
! 50: break;
! 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");
! 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_A) && (NID(b) <= N_A ) )
! 85: (*addnumt[MIN(NID(a),NID(b))])(a,b,c);
! 86: else {
! 87: if ( NID(a) != N_B ) a = tobf(a,0);
! 88: if ( NID(b) != N_B ) b = tobf(b,0);
! 89: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
! 90: mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 91: MPFRTOBF(r,d);
! 92: *c = (Num)d;
! 93: }
! 94: }
! 95:
! 96: void subbf(Num a,Num b,Num *c)
! 97: {
! 98: mpfr_t r;
! 99: BF d;
! 100:
! 101: if ( !a )
! 102: (*chsgnnumt[NID(b)])(b,c);
! 103: else if ( !b )
! 104: *c = a;
! 105: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
! 106: (*subnumt[MIN(NID(a),NID(b))])(a,b,c);
! 107: else {
! 108: if ( NID(a) != N_B ) a = tobf(a,0);
! 109: if ( NID(b) != N_B ) b = tobf(b,0);
! 110: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
! 111: mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 112: MPFRTOBF(r,d);
! 113: *c = (Num)d;
! 114: }
! 115: }
! 116:
! 117: void mulbf(Num a,Num b,Num *c)
! 118: {
! 119: mpfr_t r;
! 120: BF d;
! 121: int prec;
! 122:
! 123: if ( !a || !b )
! 124: *c = 0;
! 125: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
! 126: (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);
! 127: else {
! 128: if ( NID(a) != N_B ) a = tobf(a,0);
! 129: if ( NID(b) != N_B ) b = tobf(b,0);
! 130: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
! 131: mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 132: MPFRTOBF(r,d);
! 133: *c = (Num)d;
! 134: }
! 135: }
! 136:
! 137: void divbf(Num a,Num b,Num *c)
! 138: {
! 139: mpfr_t r;
! 140: BF d;
! 141:
! 142: if ( !b )
! 143: error("divbf : division by 0");
! 144: else if ( !a )
! 145: *c = 0;
! 146: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
! 147: (*divnumt[MIN(NID(a),NID(b))])(a,b,c);
! 148: else {
! 149: if ( NID(a) != N_B ) a = tobf(a,0);
! 150: if ( NID(b) != N_B ) b = tobf(b,0);
! 151: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
! 152: mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 153: MPFRTOBF(r,d);
! 154: *c = (Num)d;
! 155: }
! 156: }
! 157:
! 158: void pwrbf(Num a,Num b,Num *c)
! 159: {
! 160: mpfr_t r;
! 161: BF d;
! 162:
! 163: if ( !b )
! 164: *c = (Num)ONE;
! 165: else if ( !a )
! 166: *c = 0;
! 167: else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
! 168: (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c);
! 169: else {
! 170: if ( NID(a) != N_B ) a = tobf(a,0);
! 171: if ( NID(b) != N_B ) b = tobf(b,0);
! 172: mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
! 173: mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
! 174: MPFRTOBF(r,d);
! 175: *c = (Num)d;
! 176: }
! 177: }
! 178:
! 179: void chsgnbf(Num a,Num *c)
! 180: {
! 181: mpfr_t r;
! 182: BF d;
! 183:
! 184: if ( !a )
! 185: *c = 0;
! 186: else if ( NID(a) <= N_A )
! 187: (*chsgnnumt[NID(a)])(a,c);
! 188: else {
! 189: mpfr_init2(r,BFPREC(a));
! 190: mpfr_neg(r,((BF)a)->body,MPFR_RNDN);
! 191: MPFRTOBF(r,d);
! 192: *c = (Num)d;
! 193: }
! 194: }
! 195:
! 196: int cmpbf(Num a,Num b)
! 197: {
! 198: if ( !a ) {
! 199: if ( !b || (NID(b)<=N_A) )
! 200: return (*cmpnumt[NID(b)])(a,b);
! 201: else
! 202: return -mpfr_sgn(((BF)a)->body);
! 203: } else if ( !b ) {
! 204: if ( !a || (NID(a)<=N_A) )
! 205: return (*cmpnumt[NID(a)])(a,b);
! 206: else
! 207: return mpfr_sgn(((BF)a)->body);
! 208: } else {
! 209: if ( NID(a) != N_B ) a = tobf(a,0);
! 210: if ( NID(b) != N_B ) b = tobf(b,0);
! 211: return mpfr_cmp(((BF)a)->body,((BF)b)->body);
! 212: }
1.1 noro 213: }
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