[BACK]Return to bf.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

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>